RMarkdown history:

  1. Feb 2021 started github repo (https://github.com/shu251/midcayman-rise-microeuk)

  2. continue to commit changes and as you include knitr and commit changes, those will be pushed to the git repo

  3. eventually set up the gitpages under settings

  4. https://resources.github.com/whitepapers/github-and-rstudio/

1 Mid-Cayman Rise grazing

Summary of cruise operations

Below is the compilation of all microscopy data

1.1 Set up working environment

GitHub repo hosts all raw data to reproduce analysis and results. See input-data or clone the repo.

library(tidyverse); library(cowplot); library(broom)

2 Protistan cell concentration

Import eukaryotic cell count data from grazing experiments. In this section, we will calculate cells per ml from raw counts (Field of view, etc.) and use to estimate protist cell concentration. These will be used below in grazing experiment calculations.

counts <- read.delim("input-data/euk-counts-compiled.txt", 
                     blank.lines.skip = FALSE,
                     na.strings = c("", "NA"),
                     stringsAsFactors = FALSE) # Import
counts[is.na(counts)] <- 0 # Change blanks to zeroes

Raw data table collected during microscopy count process. Below code reviews the structure of this raw data and updates column headers to be more ‘R’ friendly.

colnames(counts) <- c("DATE", "SAMPLE", "EXPID", "VOL", "MAG", "FOV", "nanoNoFLP", "microNoFLP", "nanoFLP", "microFLP", "NOTES", "DateCompiled"); colnames(counts)
##  [1] "DATE"         "SAMPLE"       "EXPID"        "VOL"          "MAG"         
##  [6] "FOV"          "nanoNoFLP"    "microNoFLP"   "nanoFLP"      "microFLP"    
## [11] "NOTES"        "DateCompiled"
# Column 9 and 10 list the total number of FLP inside each, so # of occurences

# Notes if count was not countable
# unique(counts$NOTES)
head(counts)
##     DATE   SAMPLE   EXPID VOL MAG FOV nanoNoFLP microNoFLP nanoFLP microFLP
## 1 4/2/20 VD-Plume T0-Rep1 150 100   1         0          0       0        0
## 2 4/2/20 VD-Plume T0-Rep1 150 100   2         0          0       0        0
## 3 4/2/20 VD-Plume T0-Rep1 150 100   3         3          0       0        0
## 4 4/2/20 VD-Plume T0-Rep1 150 100   4         1          0       4        0
## 5 4/2/20 VD-Plume T0-Rep1 150 100   5         0          1       0        0
## 6 4/2/20 VD-Plume T0-Rep1 150 100   6         0          0       0        0
##             NOTES DateCompiled
## 1               0       6/7/20
## 2               0       6/7/20
## 3               0       6/7/20
## 4               0       6/7/20
## 5 elongated cell?       6/7/20
## 6               0       6/7/20

To count occurence and number of FLP ingested by eukaryotic cells, the number of FLPs ingested was tallied and comma separated for multiple eukaryotic cells with FLP. These values need to separated and counted as 1 eukaryotic cell each, but retain the number of FLP per cell.

counts_occur <- counts %>%
  # remove incomplete
  filter(NOTES != "Not countable") %>% 
  # Count number of euk cells observed with FLPs (ex. if "1,2", 'occur' will = 2)
  mutate(nanoFLP_occur = as.numeric(str_count(nanoFLP, "[1-9]\\d*")), 
         microFLP_occur = as.numeric(str_count(microFLP, "[1-9]\\d*")),
         # Add number of euk cells with FLPs to those without for total number of euk cells
         nanoTOTAL = as.numeric(nanoNoFLP) + nanoFLP_occur, 
         microTOTAL = as.numeric(microNoFLP) + microFLP_occur,
         euksTOTAL = nanoTOTAL + microTOTAL) %>%
      data.frame

# unique(counts$nanoFLP)
# unique(counts_occur$nanoFLP_occur)

2.1 Calculate cells per ml (euk)

Input data are the raw microscopy counts by FOV. Code below calculations cells/ml based on these values. Additionaly, variance and standard deviation are also calculated. Eukaryotic cells were also classified by size, where micro equates to >20um and nano is <20um. All counts were done at 100x magnification, confirm this:

unique(counts_occur$MAG)
## [1] 100
# head(counts_occur)

Calculate cell concentration (cells/ml).

counts_cellsml_all <- counts_occur %>%
  group_by(SAMPLE, EXPID, VOL) %>% #Calculate averages by sample
  summarise(totalFOV = n(), # Count total FOV counted
            nanoAvg = sum(nanoTOTAL)/totalFOV, #Average per FOV
            nanoVar = var(nanoTOTAL), #Variance
            nanoSd = (2*(sqrt(nanoVar))), #Standard deviation
            microAvg = sum(microTOTAL)/totalFOV, ## Repeat for microeuks
            microVar = var(microTOTAL), 
            microSd = (2*(sqrt(microVar))), 
            euksAvg = sum(euksTOTAL)/totalFOV, ## Repeat for total cell count
            euksVar = var(euksTOTAL), 
            euksSd = (2*(sqrt(euksVar))), 
            .groups = 'drop_last') %>%
  # Calculate cells/ml based on magnification (at x100, 0.01 is vol of grid), volume filtered (VOL), dilution factor (0.9), and area of counting grid (for Huber lab scope, it is 283.385):
  mutate(nanoCONC = ((nanoAvg * 283.385)/(VOL * 0.01 * 0.9)),
         microCONC = ((microAvg * 283.385)/(VOL * 0.01 * 0.9)),
         eukCONC = ((euksAvg * 283.385)/(VOL * 0.01 * 0.9))
         ) %>%
  # left_join(expmeta) %>%
  separate(SAMPLE, c("Site", "Name"), sep = "-", remove = FALSE) %>%
  separate(EXPID, c("TimePoint", "Replicate"), sep = "-", remove = FALSE) %>%
  data.frame
head(counts_cellsml_all)
##                SAMPLE    Site        Name    EXPID TimePoint Replicate VOL
## 1 Piccard-LotsOShrimp Piccard LotsOShrimp  T0-Rep3        T0      Rep3  50
## 2 Piccard-LotsOShrimp Piccard LotsOShrimp T15-Rep3       T15      Rep3  50
## 3 Piccard-LotsOShrimp Piccard LotsOShrimp T20-Rep3       T20      Rep3  50
## 4 Piccard-LotsOShrimp Piccard LotsOShrimp T40-Rep3       T40      Rep3  40
## 5       Piccard-Plume Piccard       Plume  T0-Rep1        T0      Rep1 150
## 6       Piccard-Plume Piccard       Plume  T0-Rep2        T0      Rep2 150
##   totalFOV   nanoAvg   nanoVar    nanoSd   microAvg   microVar   microSd
## 1       30 0.3666667 0.3781609 1.2298958 0.00000000 0.00000000 0.0000000
## 2       30 0.3000000 0.2172414 0.9321832 0.06666667 0.06436782 0.5074163
## 3       30 0.3333333 0.2988506 1.0933445 0.00000000 0.00000000 0.0000000
## 4       30 0.1666667 0.1436782 0.7580980 0.00000000 0.00000000 0.0000000
## 5       30 0.3000000 0.2862069 1.0699662 0.03333333 0.03333333 0.3651484
## 6       30 0.2000000 0.1655172 0.8136762 0.06666667 0.06436782 0.5074163
##     euksAvg   euksVar   euksSd  nanoCONC microCONC   eukCONC
## 1 0.3666667 0.3781609 1.229896 230.90630   0.00000 230.90630
## 2 0.3666667 0.2402299 0.980265 188.92333  41.98296 230.90630
## 3 0.3333333 0.2988506 1.093345 209.91481   0.00000 209.91481
## 4 0.1666667 0.1436782 0.758098 131.19676   0.00000 131.19676
## 5 0.3333333 0.3678161 1.212957  62.97444   6.99716  69.97160
## 6 0.2666667 0.2712644 1.041661  41.98296  13.99432  55.97728

Replicates belong to the same experiment for either Bag or IGT incubation. Below, modify these names and label new column with bag or igt. And create an average across replicates.

Average cells/ml across replicates, pivot to long format

counts_cellsml_avg <- counts_cellsml_all %>%
  select(Site, Name, TimePoint, Replicate, nanoCONC, microCONC, eukCONC) %>%
  mutate(EXP_TYPE = case_when(
    grepl("IGT", Replicate) ~ "IGT",
    grepl("Rep", Replicate) ~ "Bag"
  )) %>%
  mutate(IGT_REP = case_when(
    EXP_TYPE == "IGT" ~ Replicate,
    EXP_TYPE == "Bag" ~ "Bag")) %>%
  select(-Replicate) %>%
  pivot_longer(cols = ends_with("CONC"), names_to = "VARIABLE", values_to = "CONCENTRATION") %>%
  group_by(Site, Name, TimePoint, EXP_TYPE, IGT_REP, VARIABLE) %>%
  # Calculate mean, variance, SD, min, and max
  summarise(MEAN = mean(CONCENTRATION),
            VAR = var(CONCENTRATION),
            SD = sd(CONCENTRATION),
            SEM =(sd(CONCENTRATION)/sqrt(length(CONCENTRATION))),
            MIN = min(CONCENTRATION),
            MAX = max(CONCENTRATION),
            .groups = 'drop_last') %>%
  data.frame
head(counts_cellsml_avg)
##      Site        Name TimePoint EXP_TYPE IGT_REP  VARIABLE      MEAN VAR SD SEM
## 1 Piccard LotsOShrimp        T0      Bag     Bag   eukCONC 230.90630  NA NA  NA
## 2 Piccard LotsOShrimp        T0      Bag     Bag microCONC   0.00000  NA NA  NA
## 3 Piccard LotsOShrimp        T0      Bag     Bag  nanoCONC 230.90630  NA NA  NA
## 4 Piccard LotsOShrimp       T15      Bag     Bag   eukCONC 230.90630  NA NA  NA
## 5 Piccard LotsOShrimp       T15      Bag     Bag microCONC  41.98296  NA NA  NA
## 6 Piccard LotsOShrimp       T15      Bag     Bag  nanoCONC 188.92333  NA NA  NA
##         MIN       MAX
## 1 230.90630 230.90630
## 2   0.00000   0.00000
## 3 230.90630 230.90630
## 4 230.90630 230.90630
## 5  41.98296  41.98296
## 6 188.92333 188.92333
# View(counts_cellsml_avg)
# NOTES on calculations:
# VAR = takes the sum of the squares of each value's deviation from the mean and divides by the number of such values minus one. This differs from the calculation of variance across an entire population in that the latter divides by the size of the dataset without subtracting one.
# SD = standard deviation of all values
# SEM = standard deviation of sampling distribution; standard deviation divided by the square root of the sample size.
# Euk cell counts, all and averaged
# head(counts_cellsml_all)
# head(counts_cellsml_avg)

Parse experiment type information

# Convert to long format and add column that reports IGT vs bag experiment
plot_euk_conc <- counts_cellsml_all %>%
  select(Site, Name, TimePoint, Replicate, ends_with("CONC")) %>%
  mutate(EXP_TYPE = case_when(
    grepl("IGT", Replicate) ~ "IGT",
    grepl("Rep", Replicate) ~ "Bag"
  )) %>%
  pivot_longer(cols = ends_with("CONC"), names_to = "VARIABLE", values_to = "CONCENTRATION") %>%
  data.frame
# head(plot_euk_conc)
unique(plot_euk_conc$Name)
## [1] "LotsOShrimp"    "Plume"          "Shrimpocalypse" "BSW"           
## [5] "MustardStand"   "OMT"            "Rav2"           "ShrimpHole"    
## [9] "X18"
vent_ids <- c("BSW","Plume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole")
vent_fullname <- c("Background","Plume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole")
site_ids <- c("VD", "Piccard")
site_fullname <- c("Von Damm", "Piccard")
plot_euk_conc$SiteOrder <- factor(plot_euk_conc$Site, levels = site_ids, labels = site_fullname)
plot_euk_conc$NameOrder <- factor(plot_euk_conc$Name, levels = vent_ids, labels = vent_fullname)

2.2 Plot eukaryote cells per ml

Box plot to report all eukaryote cell counts.

# Code for plot
conc_boxplot <- ggplot(plot_euk_conc, aes(x = NameOrder, 
                                          y = CONCENTRATION, 
                                          group = NameOrder,
                                          fill = VARIABLE,
                                          shape = EXP_TYPE)) +
    geom_boxplot() + 
    # Do not color by time point
    geom_jitter(color = "black", size = 2, aes(fill = VARIABLE,
                                          shape = EXP_TYPE)) +
    scale_shape_manual(values = c(21,24)) +
    scale_fill_manual(values = c("#e7298a", "#fcbba1", "#c6dbef")) +
    coord_flip() +
    scale_y_log10() +
    # scale_y_log10(limits = c(10,1000), expand = c(0, 0)) +
    facet_grid(SiteOrder ~ EXP_TYPE, space = "free", scale = "free") +
    theme_bw() + 
  theme(axis.text.x = element_text(angle = 0, h = 1, vjust = 1),
        strip.background = element_blank(),
        legend.position = "right",
        legend.title = element_blank()) +
    labs(x = "", y = bquote("Eukaryote cells "~mL^-1),
         title = "Distribution of all eukaryotic cell counts")

Eukaryote cell concentration (cells/ml) are lower in the background and plume samples compared to vent sites. ~300 cells/ml in background and plume compared to ~1000 cells per ml at the vent sites. These values are also consistent between each vent site (Von Damm and Piccard) and between Bag and IGT samples.

Boxplot represents the median (line in box) and the 1st and 3rd quartiles in the lower and upper hinges, respectively (25th and 75th percentiles). Black data points are outliers from the boxplot. Upper and lower whiskers represent the 1.5 * interquartile ranges. Pink data points are the values contributing to the boxplot (individial counts across replicates and time points.)

conc_boxplot
## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 39 rows containing non-finite values (stat_boxplot).
## Warning: Removed 39 rows containing missing values (geom_point).

eukCONC is the sum of micro and nano. Because there was a discrepency between the micro and nano cell counts, we plan to combine for most of the analysis. Here we show that the cell concentration across replicate samples was similar throughout experiments. And that the bag versus IGT experiment results were within range of one another.

Subset to plot from T0 only.

vent_ids <- c("BSW","Plume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole")
vent_fullname <- c("Background","Plume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole")

plot_euk_format <- plot_euk_conc %>%
  filter(TimePoint == "T0" & (VARIABLE == "eukCONC")) %>%
  group_by(SiteOrder, NameOrder, TimePoint, EXP_TYPE, VARIABLE) %>%
  summarise(avg_conc = mean(CONCENTRATION),
            SEM_conc = (sd(CONCENTRATION)/sqrt(length(CONCENTRATION))),
            .groups = "rowwise") %>%
  unite(EXPERIMENT, SiteOrder, NameOrder, EXP_TYPE, remove = FALSE) %>%
  data.frame

# Factor
plot_euk_format$Site_Order <- factor(plot_euk_format$SiteOrder, levels = site_fullname, labels = site_fullname)

# View(plot_euk_format)
euk_plot <- ggplot(plot_euk_format, aes(x = NameOrder, y = avg_conc, fill = Site_Order)) +
  geom_errorbar(aes(ymax = (avg_conc + SEM_conc), ymin = (avg_conc - SEM_conc)), width = 0.2) +
  geom_point(aes(fill = Site_Order), color = "black", stat = "identity", size = 3, shape = 23) +
  facet_grid(.~ Site_Order, space = "free", scales = "free") +
  scale_fill_manual(values = c("#1c9099", "#de2d26")) +
  theme_minimal() +
    theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
           panel.background = element_blank(), 
           axis.line = element_line(colour = "black"), 
           axis.text.x = element_text(color="black", size = 12, 
                                      angle = 45, hjust = 1, vjust = 1), 
           axis.text.y = element_text(color="black", size = 12),
           axis.title =element_text(color="black", size = 12),
           axis.ticks = element_line(),
           strip.text =element_blank(), legend.title = element_blank()) +
  labs(x = "", y = bquote("Eukaryote cells "~mL^-1),
       title = "")
euk_plot

Summary of findings from eukaryote cell concentration.

# head(plot_euk_format)

plot_euk_format %>% 
  type.convert(as.is = TRUE) %>%
  filter(VARIABLE == "eukCONC") %>% 
  mutate(SAMPLE_TYPE = case_when(
    NameOrder == "Background" ~ "Background",
    NameOrder == "Plume" ~ "Plume",
    TRUE ~ SiteOrder
  )) %>%
  group_by(SAMPLE_TYPE) %>% 
  summarise(MEAN_cellml = format(mean(avg_conc), scientific = T),
           min_cellml = format(min(avg_conc), scientific = T),
           max_cellml = format(max(avg_conc), scientific = T),
           num = n())
## # A tibble: 4 × 5
##   SAMPLE_TYPE MEAN_cellml  min_cellml   max_cellml     num
##   <chr>       <chr>        <chr>        <chr>        <int>
## 1 Background  9.183773e+01 9.183773e+01 9.183773e+01     1
## 2 Piccard     3.801791e+02 2.309063e+02 4.548154e+02     3
## 3 Plume       1.185379e+02 7.930115e+01 1.577747e+02     2
## 4 Von Damm    4.105001e+02 2.597696e+02 6.20998e+02      6

3 Prokaryote cell concentration

Import counts from DAPI counts of bacteria and archaea.

prok <- read.delim("input-data/prokINSITU-counts-compiled.txt")
# head(prok)
# unique(prok$CELLML)
insitu_proks <- prok %>% 
  filter(CELLML != "not countable") %>% 
  separate(SAMPLE, c("Site", "Name"), sep = "-", remove = FALSE) %>% 
  group_by(SAMPLE, Site, Name) %>% 
  summarise(MEAN = mean(as.numeric(CELLML)),
            SD = sd(CELLML),
            SEM = (sd(CELLML)/sqrt(length(CELLML))),
            .groups = "rowwise") %>% 
  data.frame
# head(insitu_proks)
# View(insitu_proks)

3.1 Microscopy cells per ml (prok)

Plot in situ prokaryote values. Average across dives where applicable.

insitu_proks$Name_order <- factor(insitu_proks$Name, levels = c("BSW", "Plume", "Quakeplume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole", "HotChimlet1", "ShrimpGulley", "SouthofHotChimlet", "SouthofLungSnack", "ArrowLoop", "Bartizan", "Rav1"), labels = c("Background","Plume", "Quakeplume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole", "Hot Chimlet #1", "Shrimp Gulley", "South of Hot Chimlet", "South of LungSnack", "Arrow Loop", "Bartizan", "Ravelin #1"))

site_ids <- c("VD", "Piccard")
site_fullname <- c("Von Damm", "Piccard")

insitu_proks$Site_order <- factor(insitu_proks$Site, levels = site_ids, labels = site_fullname)

# unique(insitu_proks$Name)

prok_plot <- ggplot(insitu_proks, aes(x = Name_order, y = MEAN)) +
  geom_errorbar(aes(ymax = (MEAN + SEM), ymin = (MEAN - SEM)), width = 0.2) +
  geom_point(stat = "identity", shape = 23, aes(fill = Site), size = 3) +
  facet_grid(.~ Site_order, space = "free", scales = "free") +
  scale_fill_manual(values = c("#de2d26", "#1c9099")) +
  labs(y = bquote("Prokaryote cells "~mL^-1), x = "", title = "") +
  scale_y_log10() +
  theme_minimal() +
    theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
           panel.background = element_blank(), 
           axis.line = element_line(colour = "black"), 
           axis.text.x = element_text(color="black", size = 12, 
                                      angle = 45, hjust = 1, vjust = 1), 
           axis.text.y = element_text(color="black", size = 12),
           axis.title =element_text(color="black", size = 12),
           axis.ticks = element_line(),
           strip.text =element_blank(), legend.title = element_blank())
prok_plot

Summary of findings from prokaryote cell counts.

insitu_proks %>% 
  filter(Name != "Quakeplume") %>% 
  mutate(SAMPLE_TYPE = case_when(
    Name == "BSW" ~ "Background",
    Name == "Plume" ~ "Plume",
    TRUE ~ Site
  )) %>% 
  group_by(SAMPLE_TYPE) %>% 
  summarise(MEAN_cellml = format(mean(MEAN), scientific = T),
           min_cellml = format(min(MEAN), scientific = T),
           max_cellml = format(max(MEAN), scientific = T),
           num = n())
## # A tibble: 4 × 5
##   SAMPLE_TYPE MEAN_cellml  min_cellml   max_cellml     num
##   <chr>       <chr>        <chr>        <chr>        <int>
## 1 Background  2.487491e+04 1.186019e+04 3.788962e+04     2
## 2 Piccard     1.097126e+05 5.387814e+04 2.385857e+05     6
## 3 Plume       3.395372e+04 1.647831e+04 5.142913e+04     2
## 4 VD          4.090657e+04 8.816422e+03 1.114298e+05     6

3.2 Comparison with previous MCR

Compare in situ prokaryote cell counts from 2020 to previous years

prok_prev <- read.csv("input-data/cellcount_previousyr.csv")
# View(prok_prev)
# head(prok_prev)
prok_prev_formatted <- prok_prev %>% 
  mutate(VENTSITE = case_when(
    grepl("Piccard", Site) ~ "Piccard",
    grepl("Von Damm", Site) ~ "VD"
  )) %>% 
  filter(!is.na(YEAR)) %>% #QC of 
  # filter(cells_ml != "NC") %>% 
  # filter(cells_ml != "") %>% 
  # filter(cells_ml != "no data") %>% 
  type.convert(as.is = TRUE, numerals = "no.loss") %>%
  select(YEAR, VENTSITE, NAME = Name, REP=Replicate, CELLML = cells_ml, ORIGSAMPLE = Orig_vent_site_ID, ID_number, Origin)

# str(prok_prev_formatted)
# View(prok_prev_formatted)
# unique(prok_prev_formatted$VENTSITE)
# Re-import 2020
prok <- read.delim("input-data/prokINSITU-counts-compiled.txt")
# View(prok)
proks_allyrs <- prok %>% 
  separate(SAMPLE, c("VENTSITE", "NAME"), sep = "-", remove = FALSE) %>% 
  mutate(YEAR = 2020) %>%
  select(YEAR, VENTSITE, NAME, REP, CELLML, ORIGSAMPLE = BAC) %>% 
  bind_rows(prok_prev_formatted %>% select(-ID_number, -Origin)) %>% 
  type.convert(as.is = TRUE) %>%
  # Remove not countable or not data samples:
  filter(CELLML != "NC") %>%
  filter(CELLML != "") %>%
  filter(CELLML != "no data") %>%
  filter(CELLML != "not countable") %>% 
  data.frame

# View(proks_allyrs)
# View(as.data.frame(unique(proks_allyrs$NAME)))

vent_order <- c("BSW","Plume","Quakeplume","NearsummitBeebee","MainOrifice","NearMainOrifice","Rav1","HotChimlet1","HotChimlet","SouthofHotChimlet","NearHotChimlet","HotCracks1","HotCracks2","ShrimpHole","ShrimpHole(X18)","X18","X19","SouthofLungSnack","TwinPeaks","OMT","WhiteCastle","GingerCastle","ArrowLoop","Bartizan","LotsOShrimp","MustardStand","ShrimpButtery","ShrimpCanyon","ShrimpGulley","Shrimpocalypse","ShrimpVegas")
vent_names <- c("Background","Plume","Quakeplume","Near summit Beebee Vents Mound","Main Orifice","Near Main Orifice","Ravelin #1","Hot Chimlet #1","Hot Chimlet","South of Hot Chimlet","Near Hot Chimlet","Hot Cracks #1","Hot Cracks #2","Shrimp Hole","Shrimp Hole (X-18)","X-18","X-19","South of Lung Snack","Twin Peaks","Old Man Tree","White Castle","Ginger Castle","Arrow Loop","Bartizan","Lots O Shrimp","Mustard Stand","Shrimp Buttery","Shrimp Canyon","Shrimp Gulley","Shrimpocalypse","Shrimp Vegas")
proks_allyrs$NAME_ORDER <- factor(proks_allyrs$NAME, levels = vent_order, labels = vent_names)
proks_allyrs$VENTSITE_ORDER <- factor(proks_allyrs$VENTSITE, levels = c("Piccard", "VD"), labels = c("Piccard", "Von Damm"))
# pdf("compare-across-yr-cellcount-04052021.pdf", h = 8, w = 7)
ggplot(proks_allyrs, aes(x = NAME_ORDER, y = as.numeric(CELLML), fill = factor(YEAR), shape = VENTSITE_ORDER)) +
  geom_point(stat = "identity", aes(fill = factor(YEAR)), size = 3) +
  scale_shape_manual(values = c(21,23)) +
  coord_flip() +
  facet_grid(VENTSITE_ORDER ~ ., space = "free", scales = "free") +
  scale_y_log10() +
  scale_fill_manual(values = c("#1c9099", "#ffeda0", "#fc4e2a")) +
  theme_linedraw() +
  theme(axis.text = element_text(color = "black", size = 10),
        strip.background = element_blank(),
        strip.text.y = element_text(color = "black", size = 11, hjust = 0, vjust = 1),
        legend.title = element_blank(),
        legend.position = "bottom",
        panel.grid.minor = element_blank(),
        panel.grid.major = element_line(color = "grey")) +
  labs(y = bquote("Cells "~mL^-1), x = "") +
  guides(fill=guide_legend(override.aes=list(shape=22)))

# dev.off()

Generate complete table with all cell count data

# head(prok)
proks_allyrs_OUTPUT <- prok %>%
  separate(SAMPLE, c("VENTSITE", "NAME"), sep = "-", remove = FALSE) %>% 
  mutate(YEAR = 2020,
         Origin = "") %>%
  mutate(VENTSITE = case_when(
    VENTSITE == "VD" ~ "Von Damm",
    TRUE ~ VENTSITE
  )) %>% 
  separate(BAC, c("ORIGSAMPLE", "ID_number"), sep = " ") %>% 
  select(YEAR, VENTSITE, NAME, REP, CELLML, ORIGSAMPLE, ID_number, Origin) %>%
  bind_rows(prok_prev_formatted) %>%
  data.frame
## Warning: Expected 2 pieces. Additional pieces discarded in 1 rows [27].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 12 rows [17, 18,
## 19, 20, 21, 33, 34, 35, 36, 39, 40, 41].
# View(proks_allyrs_OUTPUT)
# write_delim(proks_allyrs_OUTPUT, path = "input-data/cell-counts-MCR-multiyr.txt", delim = "\t")
# head(proks_allyrs_OUTPUT)
# unique(proks_allyrs_OUTPUT$NAME)
# View(proks_allyrs_OUTPUT %>% 
#        filter(!is.na(as.numeric(CELLML))) %>% 
#        mutate(TYPE = case_when(
#          NAME == "BSW" ~ "Background",
#          NAME == "Plume" ~ "Plume",
#          TRUE ~ "Vent"
#        )) %>% 
#        mutate(VENTSITE = case_when(
#          VENTSITE == "VD" ~ "Von Damm",
#          TRUE ~ VENTSITE
#        )) %>% 
#        group_by(VENTSITE, YEAR, TYPE) %>% 
#         summarise(MEAN = mean(as.numeric(CELLML))) %>% 
#        pivot_wider(names_from = YEAR, values_from = MEAN, values_fill = NA))
# ?pivot_wider

4 Determine grazing effect

Calculate FLP per eukaryotic cell over time. Goal is to make these calculations and then determine best fit line. Slope of best fit line is the grazing rate. Need to take into account euk cells with FLPs and then the euk cells withOUT FLPs, these will be zeroes to take into account for FLPs/euk averages.

Retain previously generated dataframes:

# save(counts_cellsml_all, counts_cellsml_avg, counts_occur, file = "MCR-cellcount-dfs")
load("MCR-cellcount-dfs", verbose = TRUE)
## Loading objects:
##   counts_cellsml_all
##   counts_cellsml_avg
##   counts_occur
# Euk cell counts, all and averaged
head(counts_cellsml_all)
##                SAMPLE    Site        Name    EXPID TimePoint Replicate VOL
## 1 Piccard-LotsOShrimp Piccard LotsOShrimp  T0-Rep3        T0      Rep3  50
## 2 Piccard-LotsOShrimp Piccard LotsOShrimp T15-Rep3       T15      Rep3  50
## 3 Piccard-LotsOShrimp Piccard LotsOShrimp T20-Rep3       T20      Rep3  50
## 4 Piccard-LotsOShrimp Piccard LotsOShrimp T40-Rep3       T40      Rep3  40
## 5       Piccard-Plume Piccard       Plume  T0-Rep1        T0      Rep1 150
## 6       Piccard-Plume Piccard       Plume  T0-Rep2        T0      Rep2 150
##   totalFOV   nanoAvg   nanoVar    nanoSd   microAvg   microVar   microSd
## 1       30 0.3666667 0.3781609 1.2298958 0.00000000 0.00000000 0.0000000
## 2       30 0.3000000 0.2172414 0.9321832 0.06666667 0.06436782 0.5074163
## 3       30 0.3333333 0.2988506 1.0933445 0.00000000 0.00000000 0.0000000
## 4       30 0.1666667 0.1436782 0.7580980 0.00000000 0.00000000 0.0000000
## 5       30 0.3000000 0.2862069 1.0699662 0.03333333 0.03333333 0.3651484
## 6       30 0.2000000 0.1655172 0.8136762 0.06666667 0.06436782 0.5074163
##     euksAvg   euksVar   euksSd  nanoCONC microCONC   eukCONC
## 1 0.3666667 0.3781609 1.229896 230.90630   0.00000 230.90630
## 2 0.3666667 0.2402299 0.980265 188.92333  41.98296 230.90630
## 3 0.3333333 0.2988506 1.093345 209.91481   0.00000 209.91481
## 4 0.1666667 0.1436782 0.758098 131.19676   0.00000 131.19676
## 5 0.3333333 0.3678161 1.212957  62.97444   6.99716  69.97160
## 6 0.2666667 0.2712644 1.041661  41.98296  13.99432  55.97728
head(counts_cellsml_avg)
##      Site        Name TimePoint EXP_TYPE IGT_REP  VARIABLE      MEAN VAR SD SEM
## 1 Piccard LotsOShrimp        T0      Bag     Bag   eukCONC 230.90630  NA NA  NA
## 2 Piccard LotsOShrimp        T0      Bag     Bag microCONC   0.00000  NA NA  NA
## 3 Piccard LotsOShrimp        T0      Bag     Bag  nanoCONC 230.90630  NA NA  NA
## 4 Piccard LotsOShrimp       T15      Bag     Bag   eukCONC 230.90630  NA NA  NA
## 5 Piccard LotsOShrimp       T15      Bag     Bag microCONC  41.98296  NA NA  NA
## 6 Piccard LotsOShrimp       T15      Bag     Bag  nanoCONC 188.92333  NA NA  NA
##         MIN       MAX
## 1 230.90630 230.90630
## 2   0.00000   0.00000
## 3 230.90630 230.90630
## 4 230.90630 230.90630
## 5  41.98296  41.98296
## 6 188.92333 188.92333
head(counts_occur)
##     DATE   SAMPLE   EXPID VOL MAG FOV nanoNoFLP microNoFLP nanoFLP microFLP
## 1 4/2/20 VD-Plume T0-Rep1 150 100   1         0          0       0        0
## 2 4/2/20 VD-Plume T0-Rep1 150 100   2         0          0       0        0
## 3 4/2/20 VD-Plume T0-Rep1 150 100   3         3          0       0        0
## 4 4/2/20 VD-Plume T0-Rep1 150 100   4         1          0       4        0
## 5 4/2/20 VD-Plume T0-Rep1 150 100   5         0          1       0        0
## 6 4/2/20 VD-Plume T0-Rep1 150 100   6         0          0       0        0
##             NOTES DateCompiled nanoFLP_occur microFLP_occur nanoTOTAL
## 1               0       6/7/20             0              0         0
## 2               0       6/7/20             0              0         0
## 3               0       6/7/20             0              0         3
## 4               0       6/7/20             1              0         2
## 5 elongated cell?       6/7/20             0              0         0
## 6               0       6/7/20             0              0         0
##   microTOTAL euksTOTAL
## 1          0         0
## 2          0         0
## 3          0         3
## 4          0         2
## 5          1         1
## 6          0         0

Supplementary plot showing the concentration of eukaryotic cells over time for each experiment.

# head(counts_cellsml_avg)
vent_ids <- c("BSW","Plume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole")
vent_fullname <- c("Background","Plume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole")
site_ids <- c("VD", "Piccard")
site_fullname <- c("Von Damm", "Piccard")

counts_cellsml_avg$SiteOrder <- factor(counts_cellsml_avg$Site, levels = site_ids, labels = site_fullname)
counts_cellsml_avg$NameOrder <- factor(counts_cellsml_avg$Name, levels = vent_ids, labels = vent_fullname)

# Plot trend line of euk cell count for all experiments
counts_cellsml_avg %>%
  filter(VARIABLE == "eukCONC") %>%
  unite("Experiment", NameOrder, IGT_REP, sep = "-", remove = FALSE) %>%
  ggplot(aes(x = TimePoint, y = MEAN, shape = EXP_TYPE, fill = NameOrder)) +
    geom_path(aes(group = Experiment)) +
    # geom_errorbar(aes(ymax = (MEAN + SD), ymin = (MEAN - SD)), width = 0.2) +
    geom_errorbar(aes(ymax = (MEAN + SEM), ymin = (MEAN - SEM)), width = 0.2) +
    geom_point(stat = "identity", size = 2, aes(shape = EXP_TYPE)) +
    scale_shape_manual(values = c(21, 24)) +
    scale_fill_brewer(palette = "Paired") +
    scale_y_log10() +
    facet_wrap(SiteOrder ~ EXP_TYPE, scales = "free") +
    theme_classic() + theme(strip.background = element_blank(), 
                            legend.title = element_blank(),
                            title = element_text(size = 7, face = "bold"),
                            axis.title = element_text(size = 9)) +
    labs(title = "Total euk cell counts for each experiment", y = bquote("Average eukaryote cells "~mL^-1), x = "Time point") +
  guides(fill=guide_legend(override.aes=list(shape=21)))

## Determine FLP per euk Isolate the euk cell counts with FLPs, these are comma separated and must be separated into rows. Using counts_occur data frame from above.

# Select nano and micro counts with FLPs
counts_sepflp <- counts_occur %>% 
  filter(!NOTES == "Discard") %>% 
  filter(!(NOTES == "DTAF stain prevented counts of FLP, Euks only")) %>%
  select(DATE, SAMPLE, EXPID, VOL, MAG, FOV, nanoFLP, microFLP) %>%
  # Inputs that are comma separated will be split into a new row
  separate_rows(microFLP, sep = ",", convert = TRUE) %>%
  separate_rows(nanoFLP, sep = ",", convert = TRUE) %>%
  # Replace NAs with zeroes
  replace_na(list(microFLP = 0, nanoFLP = 0)) %>% 
  data.frame

## Check, see FOV 23, separated into rows.
# View(counts_sepflp %>%
# filter(SAMPLE == "VD-Rav2" & EXPID == "T10-Rep1"))
# View(counts_occur %>%
# filter(SAMPLE == "VD-Rav2" & EXPID == "T10-Rep1"))

Subset counts that are greater than 0, so only euk cells observed to have FLPs. Create new column that calculates FLP per Euk - by dividing by 1.

counts_flp <- counts_sepflp %>%
  select(SAMPLE, EXPID, nano_size = nanoFLP, micro_size = microFLP) %>%
  pivot_longer(cols = ends_with("_size"), names_to = "SizeFrac", values_to = "num_of_FLP") %>%
  filter(num_of_FLP > 0) %>%
  separate(SAMPLE, c("Site", "Name"), sep = "-", remove = FALSE) %>%
  separate(EXPID, c("TimePoint", "Replicate"), sep = "-", remove = FALSE) %>%
  mutate(EXP_TYPE = case_when(
    grepl("IGT", Replicate) ~ "IGT",
    grepl("Rep", Replicate) ~ "Bag"
  )) %>%
  mutate(IGT_REP = case_when(
    EXP_TYPE == "IGT" ~ Replicate,
    EXP_TYPE == "Bag" ~ "Bag")) %>%
  group_by(SAMPLE, EXPID, EXP_TYPE, IGT_REP, SizeFrac) %>%
  summarise(total_FLP = sum(num_of_FLP),
            total_euks_wflp = n(),
            .groups = "rowwise") %>%
  data.frame

OUTPUT COLUMNS: (1) total_FLP = sum of FLPs found inside a euk cell (2) total_euks_wflp = number of euks counted with ingested FLP

Repeat above operation for euk cells without any FLP. Here, subset total number of observations where there was a euk cell without FLP. These need to be counted as euk cell without an FLP.

Below code repeats process and compiles with other FLP/euk cell data.

counts_flp_compiled <- counts_occur %>% 
  filter(!(NOTES == "Discard")) %>% #Discard bad counts
  filter(!(NOTES == "DTAF stain prevented counts of FLP, Euks only")) %>%
  type.convert(as.is = TRUE) %>% #modify str() for columns
  select(SAMPLE, EXPID, nano_size = nanoNoFLP, micro_size = microNoFLP) %>% #select non flp
  pivot_longer(cols = ends_with("_size"), names_to = "SizeFrac", values_to = "num_of_euks") %>%
  separate(SAMPLE, c("Site", "Name"), sep = "-", remove = FALSE) %>%
  separate(EXPID, c("TimePoint", "Replicate"), sep = "-", remove = FALSE) %>%
  mutate(EXP_TYPE = case_when(
    grepl("IGT", Replicate) ~ "IGT",
    grepl("Rep", Replicate) ~ "Bag"
  )) %>%
  mutate(IGT_REP = case_when(
    EXP_TYPE == "IGT" ~ Replicate,
    EXP_TYPE == "Bag" ~ "Bag")) %>% 
  # filter(num_of_euks > 0) %>% # Remove observed zero counts
  group_by(SAMPLE, EXPID, EXP_TYPE, IGT_REP, SizeFrac) %>%
  summarise(total_euks_noFLP = sum(num_of_euks),
            .groups = "rowwise") %>%
  # Join with FLP count information
  ## SAMPLE, EXPID, EXPTYPE, IGTREP, and SizeFrac variables should match
  left_join(counts_flp) %>% # Join with the counts of FLP per euk cell
  replace_na(list(total_FLP = 0, total_euks_wflp = 0)) %>% #Replace NAs with zero
  data.frame
## Joining, by = c("SAMPLE", "EXPID", "EXP_TYPE", "IGT_REP", "SizeFrac")
## Check data frames
# View(counts_flp_compiled)
# View(counts_occur)
# (filter(counts_flp_compiled, SAMPLE == "VD-X18", EXPID == "T40-Rep2"))
# (filter(counts_flp_compiled, SAMPLE == "Piccard-Plume", EXPID == "T10-Rep1"))
# head(counts_flp_compiled)

Extract total euk value by adding across nano and micro, then combine back with nano and micro counts.

counts_flp_compiled_all <- counts_flp_compiled %>% 
  # Exclude size fraction:
  group_by(SAMPLE, EXPID, EXP_TYPE, IGT_REP) %>%
  summarise(total_euks_noFLP = sum(total_euks_noFLP),
            total_FLP = sum(total_FLP), 
            total_euks_wflp = sum(total_euks_wflp),
            .groups = "rowwise") %>% 
  add_column(SizeFrac = "total_euks") %>% #Add SizeFrac column
  bind_rows(counts_flp_compiled) %>% # Combine back with flp compiled list
  data.frame
# head(counts_flp_compiled_all)
# filter(counts_flp_compiled_all, SAMPLE == "VD-X18", EXPID == "T40-Rep2")
# filter(counts_flp_compiled_all, SAMPLE == "Piccard-Plume", EXPID == "T10-Rep1")

All three size fractions are represented, total, micro, and nano.

One outstanding question will be that there were many more nano size grazers than micro sized. Here I am preferring to use the total value (micro + nano), but will need to defend that this is ok.

# unique(counts_flp_compiled_all$SizeFrac)
# unique(counts_flp_compiled_all$SAMPLE)
# head(counts_flp_compiled_all[1:2,])

4.1 Calculate FLP per euk cell calculation

Import metadata which has the extact minutes for each time point used.

metadata <- read.delim("input-data/flp-exp-metadata-compiled.txt")
head(metadata)
##            SAMPLE TimePoint Minutes SiteOrigin REP
## 1 VD-MustardStand        T0       0       Vent Bag
## 2 VD-MustardStand       T10      10       Vent Bag
## 3 VD-MustardStand       T15      15       Vent Bag
## 4 VD-MustardStand       T20      20       Vent Bag
## 5 VD-MustardStand       T40      40       Vent Bag
## 6   Piccard-Plume        T0       0      Plume Bag

Add metadata information to FLP data, reformat sample names, and calculate FLP per euk as the total FLP divided by the total number of euk cells counted.

counts_flp_calcs_all <- counts_flp_compiled_all %>% 
  # Add in metadata
  # IGTXb are replicate counts, include them as replicates!
  separate(EXPID, c("TimePoint", "REP"), sep = "-", remove = FALSE) %>% mutate(
    REP = ifelse(grepl("IGT5b", REP), "IGT5", REP),
    REP = ifelse(grepl("IGT4b", REP), "IGT4", REP),
    REP = ifelse(grepl("Bag", EXP_TYPE), "Bag", REP)) %>% 
  left_join(metadata, by = c("SAMPLE" = "SAMPLE", "TimePoint" = "TimePoint", "REP" = "REP")) %>% 
  separate(SAMPLE, c("Site", "Name"), sep = "-", remove = FALSE) %>%
  separate(EXPID, c("TimePoint", "Replicate_ID"), sep = "-", remove = FALSE) %>%
  ## Treat repeated IGT counts completely separate
  # group_by(SAMPLE, Site, Name, EXPID, TimePoint, Replicate_ID, EXP_TYPE, IGT_REP, SizeFrac) %>%
  ## Treat repeated IGT counts as replicates (e.g., IGT4b and IGT4 == IGT4)
  group_by(SAMPLE, Site, Name, EXPID, TimePoint, Replicate_ID, EXP_TYPE, REP, SizeFrac) %>%
  # FLPperEuk is the total FLP divided by the total number of euk cells counted
  mutate(FLPperEuk = total_FLP/(sum(total_euks_noFLP, total_euks_wflp))) %>%
  unite("Experiment", Name, REP, sep = "-", remove = FALSE) %>%
  data.frame

# filter(counts_flp_calcs_all, SAMPLE == "VD-X18", EXPID == "T40-Rep2")
# filter(counts_flp_calcs_all, SAMPLE == "Piccard-Plume", EXPID == "T10-Rep1")

View output

# head(counts_flp_calcs_all[1:2,])
unique(counts_flp_calcs_all$Experiment)
##  [1] "LotsOShrimp-Bag"     "Plume-Bag"           "Shrimpocalypse-IGT3"
##  [4] "Shrimpocalypse-Bag"  "BSW-Bag"             "MustardStand-Bag"   
##  [7] "OMT-IGT4"            "Rav2-IGT4"           "Rav2-IGT5"          
## [10] "Rav2-Bag"            "ShrimpHole-Bag"      "X18-Bag"
# View(counts_flp_calcs_all)

COLS: Timepoint, Minutes = time point label, actual incubated minutes

COLS: Replicate_ID, REP, and IGT_REP = full replicate identified for IGTs and Bags, designation of biological replicates, and designation of technical replicates for IGT experiments

4.2 Calculate linear regression to obtain slope

Use lm() function in R to calculate linear regression for each experiment. Slope equates to grazing rate. Function inputs the FLP per euk cell data, performs regression and then adds a column for slope and r-squared values.

calculate_lm <- function(df){
  regression_1 <- df %>%
  type.convert(as.is = TRUE) %>%
  ## Keep technical replicates separate for IGTs
  # group_by(SAMPLE, Site, Experiment, Name, IGT_REP, SizeFrac) %>%
  # nest(-SAMPLE, -Site, -Experiment, -Name, -IGT_REP, -SizeFrac) %>%
  ## Combine technical replicates for IGTs
  group_by(SAMPLE, Site, Experiment, Name, REP, SizeFrac) %>%
  nest(-SAMPLE, -Site, -Experiment, -Name, -REP, -SizeFrac) %>%
  mutate(lm_fit = map(data, ~lm(FLPperEuk ~ Minutes, data = .)),
         tidied = map(lm_fit, tidy)) %>% 
  unnest(tidied) %>% 
  # select(SAMPLE, Site, Experiment, Name, IGT_REP, SizeFrac, term, estimate) %>%
  select(SAMPLE, Site, Experiment, Name, REP, SizeFrac, term, estimate) %>% 
  pivot_wider(names_from = term, values_from = estimate) %>% 
  data.frame
#   colnames(regression_1) <- c("SAMPLE", "Site", "Experiment", "Name", "IGT_REP", 
# "SizeFrac", "INTERCEPT", "SLOPE")
  colnames(regression_1) <- c("SAMPLE", "Site", 
                              "Experiment", "Name", "REP",
                              "SizeFrac", "INTERCEPT", "SLOPE")
  out_regression <- df %>%
  # group_by(SAMPLE, Site, Experiment, Name, IGT_REP, SizeFrac) %>%
  # nest(-SAMPLE, -Site, -Experiment, -Name, -IGT_REP, -SizeFrac) %>%
  group_by(SAMPLE, Site, Experiment, Name, REP, SizeFrac) %>%
  nest(-SAMPLE, -Site, -Experiment, -Name, -REP, -SizeFrac) %>%
  mutate(lm_fit = map(data, ~lm(FLPperEuk ~ Minutes, data = .)),
         glanced = map(lm_fit, glance)) %>% 
  unnest(glanced) %>% 
  # select(SAMPLE, Site, Experiment, Name, IGT_REP, SizeFrac, r.squared) %>% 
  select(SAMPLE, Site, Experiment, Name, REP, SizeFrac, r.squared) %>% 
  right_join(regression_1) %>% 
  right_join(df) %>% 
  data.frame
  out_regression$SITE <- factor(out_regression$Site, levels = c("VD", "Piccard"))
  out_regression$TYPE <- factor(out_regression$EXP_TYPE, levels = c("Bag", "IGT"))
  return(out_regression)
}

Note that an error may occur when running the below function. This is due to the fact that some experiments did not have replicates.

calcs_wslope_regression <- calculate_lm(counts_flp_calcs_all)
## Warning: All elements of `...` must be named.
## Did you want `data = c(EXPID, TimePoint, Replicate_ID, EXP_TYPE, IGT_REP, total_euks_noFLP, 
##     total_FLP, total_euks_wflp, Minutes, SiteOrigin, FLPperEuk)`?
## Warning in summary.lm(x): essentially perfect fit: summary may be unreliable
## Warning: All elements of `...` must be named.
## Did you want `data = c(EXPID, TimePoint, Replicate_ID, EXP_TYPE, IGT_REP, total_euks_noFLP, 
##     total_FLP, total_euks_wflp, Minutes, SiteOrigin, FLPperEuk)`?
## Warning in summary.lm(x): essentially perfect fit: summary may be unreliable

## Warning in summary.lm(x): essentially perfect fit: summary may be unreliable
## Joining, by = c("SAMPLE", "Site", "Experiment", "Name", "REP", "SizeFrac")
## Joining, by = c("SAMPLE", "Site", "Experiment", "Name", "REP", "SizeFrac")
# head(calcs_wslope_regression[1:2,])
# View(calcs_wslope_regression)

Use below commented out to recalculate one linear regression. Above function uses the nest() capability of tidyverse. Below, one experiment is subset to check the value.

# Extract only plume-bag experiment from VD
# tmp_plume <- filter(counts_flp_calcs_all, Experiment == "Plume-Bag") %>% filter(Site == "VD") %>% filter(SizeFrac == "total_euks")
# tmp_plume # View
# Perform linear regression
# lm_out <- lm(FLPperEuk ~ Minutes, data = tmp_plume)
# Check output
# summary(lm_out)
# lm_out$coefficients #Intercept=intercept #Minutes = SLOPE
# Compare with nested function output
# filter(calcs_wslope_regression, Experiment == "Plume-Bag") %>% filter(Site == "VD") %>% filter(SizeFrac == "total_euks") %>% head
# View(calcs_wslope_regression)
# head(calcs_wslope_regression)
# unique(calcs_wslope_regression$Site)

4.3 Plot linear regression trend

calcs_wslope_regression %>% 
  filter(SizeFrac == "total_euks") %>% 
  # filter(TYPE != "IGT") %>% 
  unite(EXPERIMENT, SITE, Experiment, sep = " ", remove = FALSE) %>% 
  ggplot(aes(x = Minutes, y = FLPperEuk, fill = Site, shape = TYPE, group = Experiment)) +
  geom_abline(aes(slope = SLOPE, intercept = INTERCEPT), color = "black", linetype = "dashed", size = 1) +
  geom_point(stat = "identity", color = "black", 
             size = 2, aes(shape = TYPE, fill = Site)) +
  scale_shape_manual(values = c(21, 24)) +
  scale_fill_manual(values = c("#de2d26", "#1c9099")) +
  labs(x = "Minutes", y = bquote("FLP"~eukaryote^-1), title = "Grazing experiment regression") +
  facet_wrap(. ~ EXPERIMENT) +
  # Report r.squared
  geom_text(aes(x = 42, y = max(FLPperEuk), label = paste(round(SLOPE, 4))), 
            vjust = 1, hjust = 0, size = 3) +
  theme_bw() + 
  theme(strip.background = element_blank(),
        strip.text = element_text(color = "black", size = 7),
                     legend.title = element_blank(),
                     legend.position = "right")

Data points represent the FLP per euk cells (based on total eukaryote cells counts). Y-axis represents the duration of incubation (in minutes). The dashed purple line reprents the slope and intercept of the experiment.

4.4 Remove Tf-IGT experiments

IGT experiment results appear to have bottle effect, especially in the final time point. Additionally, due to the lack of biological replicates in the IGT experiments, technical replicates are treated as biological replicates in the regression below.

IGT_lm_woTf <- counts_flp_calcs_all %>% 
  # Select only IGT experiments with total eukaryotes, remove Tf (T3)
  filter(SizeFrac == "total_euks") %>% 
  filter(EXP_TYPE == "IGT" & !(TimePoint == "T3")) %>% 
  add_column(IGT_cor = "rm Tf") %>% 
  data.frame

# Recalculate lm(), keep replicates separate
igt_regression_noTf <- calculate_lm(IGT_lm_woTf) # Recalculate
## Warning: All elements of `...` must be named.
## Did you want `data = c(EXPID, TimePoint, Replicate_ID, EXP_TYPE, IGT_REP, total_euks_noFLP, 
##     total_FLP, total_euks_wflp, Minutes, SiteOrigin, FLPperEuk, 
##     IGT_cor)`?

## Warning: All elements of `...` must be named.
## Did you want `data = c(EXPID, TimePoint, Replicate_ID, EXP_TYPE, IGT_REP, total_euks_noFLP, 
##     total_FLP, total_euks_wflp, Minutes, SiteOrigin, FLPperEuk, 
##     IGT_cor)`?
## Joining, by = c("SAMPLE", "Site", "Experiment", "Name", "REP", "SizeFrac")
## Joining, by = c("SAMPLE", "Site", "Experiment", "Name", "REP", "SizeFrac")

Plot IGT grazing experiment with newly calculated grazing effect.

igt_regression_noTf %>% 
  # filter(SizeFrac == "total_euks") %>% 
  # filter(TYPE != "IGT") %>% 
  unite(EXPERIMENT, SITE, Experiment, sep = " ", remove = FALSE) %>% 
  ggplot(aes(x = Minutes, y = FLPperEuk, fill = Site, shape = TYPE, group = Experiment)) +
  geom_abline(aes(slope = SLOPE, intercept = INTERCEPT), color = "black", linetype = "dashed", size = 1) +
  geom_point(stat = "identity", color = "black", 
             size = 2, aes(shape = TYPE, fill = Site)) +
  scale_shape_manual(values = c(24)) +
  scale_fill_manual(values = c("#de2d26", "#1c9099")) +
  labs(x = "Minutes", y = bquote("FLP"~eukaryote^-1), title = "Grazing experiment regression") +
  facet_wrap(. ~ EXPERIMENT) +
  # Report r.squared
  geom_text(aes(x = 5, y = max(FLPperEuk), label = paste(round(SLOPE, 4))), 
            vjust = 1, hjust = 0, size = 3) +
  theme_bw() + 
  theme(strip.background = element_blank(),
        strip.text = element_text(color = "black", size = 7),
                     legend.title = element_blank(),
                     legend.position = "right")

IGT grazing regression without the Tf data point appear more consistent. Keeping that one.

4.5 Compile grazing experiment results

calcs_wslope_regression_update <- calcs_wslope_regression %>% 
  filter(TYPE != "IGT") %>% 
  bind_rows(igt_regression_noTf %>% select(-IGT_cor)) %>% 
  data.frame

# Factor
vent_ids <- c("BSW","Plume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole")
vent_fullname <- c("Background","Plume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole")
site_ids <- c("VD", "Piccard")
site_fullname <- c("Von Damm", "Piccard")
# Factor for shipboard
calcs_wslope_regression_update$SiteOrder <- factor(calcs_wslope_regression_update$Site, levels = site_ids, labels = site_fullname)
calcs_wslope_regression_update$NameOrder <- factor(calcs_wslope_regression_update$Name, levels = vent_ids, labels = vent_fullname)

# dim(calcs_wslope_regression); dim(calcs_wslope_regression_update)
head(calcs_wslope_regression_update)
##                SAMPLE    Site      Experiment        Name REP   SizeFrac
## 1 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 2 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 3 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 4 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 5       Piccard-Plume Piccard       Plume-Bag       Plume Bag total_euks
## 6       Piccard-Plume Piccard       Plume-Bag       Plume Bag total_euks
##    r.squared INTERCEPT        SLOPE    EXPID TimePoint Replicate_ID EXP_TYPE
## 1 0.66025258 0.3517002 -0.007605829  T0-Rep3        T0         Rep3      Bag
## 2 0.66025258 0.3517002 -0.007605829 T15-Rep3       T15         Rep3      Bag
## 3 0.66025258 0.3517002 -0.007605829 T20-Rep3       T20         Rep3      Bag
## 4 0.66025258 0.3517002 -0.007605829 T40-Rep3       T40         Rep3      Bag
## 5 0.02810998 0.7435088  0.005362856  T0-Rep1        T0         Rep1      Bag
## 6 0.02810998 0.7435088  0.005362856  T0-Rep2        T0         Rep2      Bag
##   IGT_REP total_euks_noFLP total_FLP total_euks_wflp Minutes SiteOrigin
## 1     Bag                9         3               2       0       Vent
## 2     Bag                8         4               3      15       Vent
## 3     Bag                9         2               1      20       Vent
## 4     Bag                5         0               0      40       Vent
## 5     Bag                6         6               4       0      Plume
## 6     Bag                5         5               3       0      Plume
##   FLPperEuk    SITE TYPE SiteOrder      NameOrder
## 1 0.2727273 Piccard  Bag   Piccard Lots 'O Shrimp
## 2 0.3636364 Piccard  Bag   Piccard Lots 'O Shrimp
## 3 0.2000000 Piccard  Bag   Piccard Lots 'O Shrimp
## 4 0.0000000 Piccard  Bag   Piccard Lots 'O Shrimp
## 5 0.6000000 Piccard  Bag   Piccard          Plume
## 6 0.6250000 Piccard  Bag   Piccard          Plume

4.6 Check FLP control experiments

All incubations had control experiments run alongside them. This was to ensure added FLP did not decrease or change in concentration over time.

bac_ctrl <- read.delim("input-data/bac-counts-compiled.txt")
# dim(bac_ctrl)
dtaf <- bac_ctrl %>% 
  separate(SampleID, c("exp", "Replicate", "TimePoint"), sep = "-", remove = FALSE) %>% 
  separate(Site, c("Site", "Name"), sep = "-", remove = FALSE) %>% 
  filter(Stain == "DTAF") %>% 
  data.frame
## Warning: Expected 2 pieces. Additional pieces discarded in 17 rows [33, 34, 35,
## 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49].
# View(bac_ctrl)
# head(dtaf)

dtaf_avg <- dtaf %>% 
  group_by(TimePoint, Stain, Site, Name) %>% 
  summarise(Avg_cellsperml = mean(Cells.ml)) %>% 
  data.frame
## `summarise()` has grouped output by 'TimePoint', 'Stain', 'Site'. You can override using the `.groups` argument.
dtaf_avg %>% 
  filter(Site != "IGT") %>% 
  ggplot(aes(x = TimePoint, y = Avg_cellsperml, fill = Name, shape = Site)) +
  geom_rect(data = filter(dtaf_avg, TimePoint == "T0", Site != "IGT"), aes(
                                           ymin = (Avg_cellsperml-(0.1*Avg_cellsperml)),
                                           ymax = (Avg_cellsperml+(0.1*Avg_cellsperml))), color = NA, alpha = 0.4, xmin = 0, xmax = 6, fill = "black") +
  geom_line(aes(group = Name)) +
  geom_point(stat = "identity", aes(shape = Site, fill = Name), size = 2) +
  # scale_fill_manual(values = c("black","#9970ab", "#5aae61")) +
  facet_wrap(Name ~ Site) +
  scale_y_log10() +
  theme_bw() + theme(strip.background = element_blank(), 
                          legend.title = element_blank(),
                     axis.text = element_text(size = 10, color = "black"),
                          title = element_text(size = 10, face = "bold"),
                          axis.title = element_text(size = 9)) +
  labs(title = "FLP counts for controls", y = bquote("Log FLP "~mL^-1), x = "Time point")

Control experiments look ok for now.

dtaf_avg %>% 
  filter(Site == "IGT") %>% 
  ggplot(aes(x = TimePoint, y = Avg_cellsperml, fill = Name, shape = Site)) +
  geom_rect(data = filter(dtaf_avg, TimePoint == "T0", Site == "IGT"), aes(
                                           ymin = (Avg_cellsperml-(0.1*Avg_cellsperml)),
                                           ymax = (Avg_cellsperml+(0.1*Avg_cellsperml))), color = NA, alpha = 0.4, xmin = 0, xmax = 6, fill = "black") +
  geom_line(aes(group = Name)) +
  geom_point(stat = "identity", aes(shape = Site, fill = Name), size = 2) +
  # scale_fill_manual(values = c("black","#9970ab", "#5aae61")) +
  facet_wrap(Name ~ Site) +
  scale_y_log10() +
  theme_bw() + theme(strip.background = element_blank(), 
                          legend.title = element_blank(),
                     axis.text = element_text(size = 10, color = "black"),
                          title = element_text(size = 10, face = "bold"),
                          axis.title = element_text(size = 9)) +
  labs(title = "FLP counts for controls", y = bquote("Log FLP "~mL^-1), x = "Time point")

4.7 Table with grazing values

head(calcs_wslope_regression_update)
##                SAMPLE    Site      Experiment        Name REP   SizeFrac
## 1 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 2 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 3 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 4 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 5       Piccard-Plume Piccard       Plume-Bag       Plume Bag total_euks
## 6       Piccard-Plume Piccard       Plume-Bag       Plume Bag total_euks
##    r.squared INTERCEPT        SLOPE    EXPID TimePoint Replicate_ID EXP_TYPE
## 1 0.66025258 0.3517002 -0.007605829  T0-Rep3        T0         Rep3      Bag
## 2 0.66025258 0.3517002 -0.007605829 T15-Rep3       T15         Rep3      Bag
## 3 0.66025258 0.3517002 -0.007605829 T20-Rep3       T20         Rep3      Bag
## 4 0.66025258 0.3517002 -0.007605829 T40-Rep3       T40         Rep3      Bag
## 5 0.02810998 0.7435088  0.005362856  T0-Rep1        T0         Rep1      Bag
## 6 0.02810998 0.7435088  0.005362856  T0-Rep2        T0         Rep2      Bag
##   IGT_REP total_euks_noFLP total_FLP total_euks_wflp Minutes SiteOrigin
## 1     Bag                9         3               2       0       Vent
## 2     Bag                8         4               3      15       Vent
## 3     Bag                9         2               1      20       Vent
## 4     Bag                5         0               0      40       Vent
## 5     Bag                6         6               4       0      Plume
## 6     Bag                5         5               3       0      Plume
##   FLPperEuk    SITE TYPE SiteOrder      NameOrder
## 1 0.2727273 Piccard  Bag   Piccard Lots 'O Shrimp
## 2 0.3636364 Piccard  Bag   Piccard Lots 'O Shrimp
## 3 0.2000000 Piccard  Bag   Piccard Lots 'O Shrimp
## 4 0.0000000 Piccard  Bag   Piccard Lots 'O Shrimp
## 5 0.6000000 Piccard  Bag   Piccard          Plume
## 6 0.6250000 Piccard  Bag   Piccard          Plume
# Generate final table
table_grazerate <- calcs_wslope_regression_update %>% 
  filter(SizeFrac == "total_euks") %>% 
  group_by(SAMPLE, SiteOrder, NameOrder, SLOPE, EXP_TYPE) %>% 
  mutate(No_of_Replicates = n_distinct(Replicate_ID)) %>%
  select(SAMPLE, Site=SiteOrder, Name=NameOrder, RATE = SLOPE, EXP_TYPE, No_of_Replicates) %>% 
  distinct() %>% 
  data.frame
# View(table_grazerate)
# write_delim(table_grazerate, path = "Table-grazingrates.txt", delim = "\t")
library(gt)
table_grazerate %>% 
  group_by(Site, EXP_TYPE) %>% 
  gt() #%>% 
SAMPLE Name RATE No_of_Replicates
Piccard - Bag
Piccard-LotsOShrimp Lots 'O Shrimp -0.007605829 1
Piccard-Plume Plume 0.005362856 3
Piccard-Shrimpocalypse Shrimpocalypse 0.015686872 1
Von Damm - Bag
VD-BSW Background 0.002958889 3
VD-MustardStand Mustard Stand -0.005445545 2
VD-Plume Plume 0.005274231 3
VD-Rav2 Ravelin #2 0.003470217 3
VD-ShrimpHole Shrimp Hole -0.001967253 2
VD-X18 X-18 0.001744429 2
Piccard - IGT
Piccard-Shrimpocalypse Shrimpocalypse 0.016794258 1
Von Damm - IGT
VD-OMT Old Man Tree 0.014510943 2
VD-Rav2 Ravelin #2 0.015395240 2
VD-Rav2 Ravelin #2 0.001603035 2
  # fmt_number(columns = 3, decimals = 4)

Add FLP conc to table

# head(table_grazerate)
dtaf_igt <- 5352.8278 # Manually insert FLP concentration for IGT experiments; this value is estimated from how IGT FLP spike-ins were calculated

table_grazerate_wflp <- bac_ctrl %>% 
  filter(FLP_t0 == "use") %>% 
  add_column(EXP_TYPE = "Bag") %>% 
  group_by(Site, EXP_TYPE) %>% 
  summarise(FLP_conc = mean(Cells.ml)) %>% 
  right_join(table_grazerate, by = c("Site" = "SAMPLE", "EXP_TYPE" = "EXP_TYPE")) %>% 
  mutate(FLP_conc = ifelse(EXP_TYPE == "IGT", dtaf_igt, FLP_conc)) %>% 
  data.frame
## `summarise()` has grouped output by 'Site'. You can override using the `.groups` argument.
# View(table_grazerate_wflp)

Negative grazing rates equal undetected (change to zero).

# head(table_grazerate)
bsw <- c("Plume", "Background")
grazerate_plot <- table_grazerate %>% 
  type.convert(as.is = TRUE) %>% 
  mutate(detected = case_when(
    RATE < 0 ~ "Not detected",
    TRUE ~ "Detected")) %>% 
  mutate(type = case_when(
    Name %in% bsw ~ Name,
    TRUE ~ paste("Vent", EXP_TYPE, sep="-")
  )) %>% 
  mutate(GRAZE_RATE = case_when(
    RATE < 0 ~ 0,
    TRUE ~ RATE
  )) %>% 
  mutate(type_site = case_when(
    Name %in% bsw ~ Name,
    TRUE ~ "Vent"
  )) %>% 
  # filter(detected == "Detected") %>% 
  data.frame
# Factoring
type_order <- c("Vent-Bag", "Vent-IGT", "Plume", "Background")
grazerate_plot$TYPE <- factor(grazerate_plot$type, levels = type_order)
vent_ids <- c("BSW","Plume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole")
vent_fullname <- c("Background","Plume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole")
site_ids <- c("VD", "Piccard")
site_fullname <- c("Von Damm", "Piccard")
vent_colors <- c("#0868ac", "#41ab5d", "#e7298a", "#c994c7", "#fc4e2a", "#fed976", "#6a51a3", "#ffeda0", "#a1d99b")
names(vent_colors) <- vent_fullname
grazerate_plot$NAME <- factor(grazerate_plot$Name, levels = vent_fullname)
# svg("", h =, w = )
grazing_min_plot <- grazerate_plot %>% 
  ggplot(aes(y = GRAZE_RATE, x = NAME, shape = EXP_TYPE, fill = Site)) +
  geom_jitter(stat = "identity", aes(shape = EXP_TYPE, fill = Site),
              color = "black", size = 3, width = 0.3) +
  scale_shape_manual(values = c(21, 24)) +
  scale_fill_manual(values = c("#de2d26", "#1c9099")) +
  facet_grid(.~Site, space = "free", scales = "free") +
  # coord_flip() +
    theme_minimal() +
    theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
           panel.background = element_blank(), 
           axis.line = element_line(colour = "black"), 
           axis.text.x = element_text(color="black", size = 12, 
                                      angle = 45, hjust = 1, vjust = 1), 
           axis.text.y = element_text(color="black", size = 12),
           axis.title =element_text(color="black", size = 12),
           axis.ticks = element_line(),
           strip.text =element_blank(), legend.title = element_blank())+
    guides(fill = guide_legend(override.aes = list(shape = c(21))),
       shape = guide_legend(override.aes = list(fill = "black"))) +
    labs(x = "", y = bquote("FLP " ~grazer^-1 ~min^-1))
# dev.off()
grazing_min_plot

5 Plot cell conc & grazing

Generate final plots that show eukaryote and prokaryote biomass and cell concentration. Save final table reporting carbon estimates.

# Subset the average in situ prok cells/ml for non-background samples
tmp <- filter(insitu_proks, Name != "BSW", Name != "Plume") %>% select(MEAN)
avg_insitu <- mean(tmp$MEAN)
# head(insitu_proks)

# Add to master table with data
table_grazerate_wflp_wprok <- insitu_proks %>% 
  select(Site = SAMPLE, Prok_conc = MEAN, Prok_sem = SEM) %>% 
  right_join(table_grazerate_wflp) %>% 
  mutate(Prok_conc = ifelse(is.na(Prok_conc), avg_insitu, Prok_conc)) %>% 
  data.frame
## Joining, by = "Site"

Table with all cell count data

# head(table_grazerate_wflp_wprok)
table_grazerate_wflp_wprok_weuk <- plot_euk_format %>% 
  select(Name = NameOrder, Site.y.y = SiteOrder, euk_conc = avg_conc, EXP_TYPE, euk_conc_sem = SEM_conc) %>% 
  right_join(table_grazerate_wflp_wprok) %>% 
  select(SITE = Site.y.y, NAME = Name, EXP = EXP_TYPE, SAMPLE = Site, RATE_min = RATE, REPs = No_of_Replicates, FLP_ml = FLP_conc, PROK_ml = Prok_conc, PROK_sem = Prok_sem, EUK_ml = euk_conc, EUK_sem = euk_conc_sem) %>% 
  mutate(RATE_min = case_when(
    RATE_min < 0 ~ 0,
    TRUE ~ RATE_min
  )) %>% 
  data.frame
## Joining, by = c("Name", "Site.y.y", "EXP_TYPE")

6 Grazing rate calculations

# head(table_grazerate_wflp_wprok_weuk)

Based on Unrein et al. 2007, we use the estimated grazing rate, in situ prok abundance, in situ euk abundance, and the concentration of FLP to make additional estimates.

SKH: add citations for the below calulcations.

table_wcalcs <- table_grazerate_wflp_wprok_weuk %>%
  # Ingestion rate per hour
  mutate(RATE_hr = (RATE_min * 60),
         RATE_day = (RATE_hr * 24), #Compare to GR?
         # FLP concentration per L
         FLP_L = (FLP_ml * 1000),
         # mL per grazer per hr
         CLEARANCE_RATE_ml = (RATE_hr/FLP_ml),
         # nL per grazer per hour
         CLEARANCE_RATE_nL = ((RATE_hr/FLP_ml)/1.00E+6), 
         # proks per grazer per hr
         SPEC_GRAZE_RATE_hr = (CLEARANCE_RATE_ml * PROK_ml), 
         # proks per grazer per day
         GRAZE_RATE_DAY = (24 * SPEC_GRAZE_RATE_hr),
         # proks per ml per hr
         GRAZING_EFFECT_hr = (SPEC_GRAZE_RATE_hr * EUK_ml),
         GRAZING_EFFECT_hr_min = (SPEC_GRAZE_RATE_hr * (EUK_ml - EUK_sem)),
         GRAZING_EFFECT_hr_max = (SPEC_GRAZE_RATE_hr * (EUK_ml + EUK_sem)),
         # cells per ml per day
         GRAZING_EFFECT_day = ((SPEC_GRAZE_RATE_hr * 24) * EUK_ml),
         # Percentage per day
         BAC_TURNOVER_PERC = 100*(GRAZING_EFFECT_day / PROK_ml),
         BAC_TURNOVER_PERC_min = 100*(GRAZING_EFFECT_day / (PROK_ml - PROK_sem)),
         BAC_TURNOVER_PERC_max = 100*(GRAZING_EFFECT_day / (PROK_ml + PROK_sem))) %>% 
  data.frame
# View(table_wcalcs)

Explanation of units for table with calculated values.

  • RATE_min & RATE_hr = Grazing rate as ‘FLPs per grazer per minute’ and per hour

  • CLEARANCE_RATE = ml or nL per grazer per hour

  • SPEC_GRAZE_RATE (Specific grazing rate) = Prokaryotes per grazer per hour

  • GRAZING EFFECT = bacteria per ml per hour

  • Bacterial turnover rate = % per day

7 Carbon biomass table construction

Unrein, F., Massana, R., Alonso-Sáez, L., and Gasol, J.M. (2007) Significant year-round effect of small mixotrophic flagellates on bacterioplankton in an oligotrophic coastal system. Limnol Oceanogr 52: 456–469.

# colnames(table_wcalcs)
bkgd <- c("Background", "Plume")
table_wcalcs %>% 
  mutate(loc_type = case_when(
    NAME %in% bkgd ~ "Background",
    TRUE ~ "Vent fluid"
  )) %>%
  # group_by(loc_type, SITE, EXP) %>% 
  select(-SAMPLE) %>% 
  gt(
    groupname_col = c("SITE", "EXP", "loc_type"),
    rowname_col = "NAME"
  ) %>% 
  cols_label(RATE_min = html("minute<sup>-1</sup>"),
             RATE_hr = html("hour<sup>-1</sup>"),
             RATE_day = html("day<sup>-1</sup>"),
             REPs = html("# of incubations"),
             FLP_ml = html("FLP ml<sup>-1</sup>"),
             PROK_ml = html("Prokaryote cells ml<sup>-1</sup>"),
             PROK_sem = html("SEM prokaryote cells ml<sup>-1</sup>"),
             EUK_ml = html("Eukaryote cells ml<sup>-1</sup>"),
             EUK_sem = html("SEM eukaryote cells ml<sup>-1</sup>"),
             FLP_L = html("FLP L<sup>-1</sup>"),
             CLEARANCE_RATE_ml = html("ml grazer<sup>-1</sup> hr<sup>-1</sup>"),
             CLEARANCE_RATE_nL = html("nl grazer<sup>-1</sup> hr<sup>-1</sup>"),
             SPEC_GRAZE_RATE_hr = html("Prokaryote grazer<sup>-1</sup> hr<sup>-1</sup>"),
             GRAZE_RATE_DAY = html("Prokaryote grazer<sup>-1</sup> day<sup>-1</sup>"),
             GRAZING_EFFECT_hr = html("Prokaryote ml<sup>-1</sup> hr<sup>-1</sup>"),
             GRAZING_EFFECT_hr_min = html("MIN"),
             GRAZING_EFFECT_hr_max = html("MAX"),
             GRAZING_EFFECT_day = html("Prokaryote ml<sup>-1</sup> day<sup>-1</sup>"),
             BAC_TURNOVER_PERC = html("Bacteria turnover % day<sup>-1</sup>"),
             BAC_TURNOVER_PERC_min = html("MIN"),
             BAC_TURNOVER_PERC_max = html("MAX")) %>% 
  tab_spanner(
    label = (html("Turnover")),
    columns = starts_with("BAC_TURNOVER")
  ) %>% 
  tab_spanner(
    label = (html("Grazing rate: prokaryote cells consumed")),
    columns = starts_with("GRAZING_EFFECT")
  ) %>% 
  tab_spanner(
    label = (html("ml grazer<sup>-1</sup> hr<sup>-1</sup>")),
    columns = c(CLEARANCE_RATE_ml, CLEARANCE_RATE_nL)
  ) %>% 
  tab_spanner(
    label = html("Specific grazing rate"),
    columns = c(SPEC_GRAZE_RATE_hr, GRAZE_RATE_DAY)
  ) %>% 
  tab_spanner(
    label = (html("FLPs grazer<sup>-1</sup>")),
    columns = c(RATE_hr, RATE_min, RATE_day),
  ) %>% 
  tab_spanner(
    label = (html("Cell counts")),
    columns = c(PROK_ml, PROK_sem, EUK_ml, EUK_sem, FLP_L, FLP_ml),
  ) %>% 
  tab_source_note(source_note = "NAs indicate values were unavailable.
                  Zero values for rates indicate no grazing pressure detected.") %>% 
  fmt_scientific(columns = everything()) %>% 
  tab_options(
    table.font.size = 12,
    table.border.top.color = "black",
    column_labels.border.bottom.color = "black",
    column_labels.border.bottom.width= px(3),
    table.width = pct(100))
# of incubations Cell counts FLPs grazer-1 ml grazer-1 hr-1 Specific grazing rate Grazing rate: prokaryote cells consumed Turnover
Prokaryote cells ml-1 SEM prokaryote cells ml-1 Eukaryote cells ml-1 SEM eukaryote cells ml-1 FLP L-1 FLP ml-1 hour-1 minute-1 day-1 ml grazer-1 hr-1 nl grazer-1 hr-1 Prokaryote grazer-1 hr-1 Prokaryote grazer-1 day-1 Prokaryote ml-1 hr-1 MIN MAX Prokaryote ml-1 day-1 Bacteria turnover % day-1 MIN MAX
Von Damm - Bag - Background
Background 3.00 3.79 × 104 8.61 × 103 9.18 × 101 2.19 × 101 4.86 × 106 4.86 × 103 1.78 × 10−1 2.96 × 10−3 4.26 3.65 × 10−5 3.65 × 10−11 1.38 3.32 × 101 1.27 × 102 9.68 × 101 1.57 × 102 3.05 × 103 8.05 1.04 × 101 6.56
Plume 3.00 1.65 × 104 2.62 × 103 1.58 × 102 6.71 × 101 3.42 × 107 3.42 × 104 3.16 × 10−1 5.27 × 10−3 7.59 9.24 × 10−6 9.24 × 10−12 1.52 × 10−1 3.65 2.40 × 101 1.38 × 101 3.42 × 101 5.77 × 102 3.50 4.16 3.02
Von Damm - Bag - Vent fluid
X-18 2.00 1.11 × 105 2.97 × 103 3.15 × 102 1.05 × 102 3.15 × 106 3.15 × 103 1.05 × 10−1 1.74 × 10−3 2.51 3.32 × 10−5 3.32 × 10−11 3.70 8.89 × 101 1.17 × 103 7.78 × 102 1.56 × 103 2.80 × 104 2.51 × 101 2.58 × 101 2.45 × 101
Ravelin #2 3.00 7.11 × 104 NA 4.09 × 102 7.35 × 101 5.19 × 107 5.19 × 104 2.08 × 10−1 3.47 × 10−3 5.00 4.01 × 10−6 4.01 × 10−12 2.85 × 10−1 6.85 1.17 × 102 9.59 × 101 1.38 × 102 2.80 × 103 3.94 NA NA
Mustard Stand 2.00 5.67 × 104 1.43 × 104 2.60 × 102 2.89 × 101 6.60 × 106 6.60 × 103 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
Shrimp Hole 2.00 4.20 × 104 3.15 × 103 3.86 × 102 7.87 1.13 × 108 1.13 × 105 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
Von Damm - IGT - Vent fluid
Old Man Tree 2.00 7.11 × 104 NA 4.72 × 102 1.22 × 102 5.35 × 106 5.35 × 103 8.71 × 10−1 1.45 × 10−2 2.09 × 101 1.63 × 10−4 1.63 × 10−10 1.16 × 101 2.78 × 102 5.47 × 103 4.05 × 103 6.88 × 103 1.31 × 105 1.84 × 102 NA NA
Ravelin #2 2.00 7.11 × 104 NA 6.21 × 102 1.23 × 102 5.35 × 106 5.35 × 103 9.24 × 10−1 1.54 × 10−2 2.22 × 101 1.73 × 10−4 1.73 × 10−10 1.23 × 101 2.95 × 102 7.62 × 103 6.11 × 103 9.14 × 103 1.83 × 105 2.57 × 102 NA NA
Ravelin #2 2.00 7.11 × 104 NA 6.21 × 102 1.23 × 102 5.35 × 106 5.35 × 103 9.62 × 10−2 1.60 × 10−3 2.31 1.80 × 10−5 1.80 × 10−11 1.28 3.07 × 101 7.94 × 102 6.36 × 102 9.51 × 102 1.91 × 104 2.68 × 101 NA NA
Piccard - Bag - Background
Plume 3.00 5.14 × 104 4.62 × 103 7.93 × 101 1.68 × 101 2.97 × 107 2.97 × 104 3.22 × 10−1 5.36 × 10−3 7.72 1.09 × 10−5 1.09 × 10−11 5.58 × 10−1 1.34 × 101 4.43 × 101 3.49 × 101 5.36 × 101 1.06 × 103 2.07 2.27 1.90
Piccard - Bag - Vent fluid
Shrimpocalypse 1.00 2.39 × 105 6.58 × 104 4.55 × 102 NA 1.70 × 107 1.70 × 104 9.41 × 10−1 1.57 × 10−2 2.26 × 101 5.54 × 10−5 5.54 × 10−11 1.32 × 101 3.17 × 102 6.01 × 103 NA NA 1.44 × 105 6.04 × 101 8.34 × 101 4.74 × 101
Lots 'O Shrimp 1.00 5.39 × 104 1.37 × 104 2.31 × 102 NA 6.17 × 105 6.17 × 102 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 NA NA 0.00 0.00 0.00 0.00
Piccard - IGT - Vent fluid
Shrimpocalypse 1.00 2.39 × 105 6.58 × 104 4.55 × 102 7.00 × 101 5.35 × 106 5.35 × 103 1.01 1.68 × 10−2 2.42 × 101 1.88 × 10−4 1.88 × 10−10 4.49 × 101 1.08 × 103 2.04 × 104 1.73 × 104 2.36 × 104 4.90 × 105 2.05 × 102 2.84 × 102 1.61 × 102
NAs indicate values were unavailable. Zero values for rates indicate no grazing pressure detected.

8 Estimates of carbon biomass and consumption

8.1 Protist carbon-measured

References for estimating biovolume Pernice, M.C., Forn, I., Gomes, A., Lara, E., Alonso-Sáez, L., Arrieta, J.M., et al. (2015) Global abundance of planktonic heterotrophic protists in the deep ocean. ISME J 9: 782–792.

# Import manual biovolume measurements
biov <- read.delim("input-data/biovol-euk-12-10-2020.txt")
# head(biov)

# Calculate volume
biov_calc <- biov %>% 
  mutate(SizeFrac = case_when(
    h >= 20 ~ "micro",
    TRUE ~ "nano")) %>% 
  mutate(Volume = ((pi/6) * (d^2) * d)) %>% # Calculate volume (um cubed) # Hillebrand et al. 1999
  mutate(pgC_cell = (0.216 * (Volume^0.939))) %>% # Calculate Cell biomass in pg C per cell # Menden-Deuer and Lessard 2000
  data.frame
# View(biov_calc)
biov_calc
##    EXP VENT_BSW      h      d SizeFrac     Volume     pgC_cell
## 1  IGT     vent 30.077 25.764    micro 8954.44130 1110.2426245
## 2  IGT     vent 89.582 10.000    micro  523.59878   77.1957618
## 3  Bag      BSW 14.595  8.036     nano  271.71800   41.6956679
## 4  Bag      BSW 12.480  8.982     nano  379.41786   57.0486460
## 5  Bag     vent  9.218  3.120     nano   15.90239    2.9015292
## 6  IGT     vent 17.255  9.986     nano  521.40274   76.8917043
## 7  Bag     vent 41.153 21.000    micro 4849.04826  624.1445904
## 8  IGT     vent 10.282  4.136     nano   37.04591    6.4194942
## 9  IGT     vent 29.776 25.852    micro 9046.50993 1120.9583343
## 10 IGT     vent 10.991  4.000     nano   33.51032    5.8424695
## 11 Bag     vent 14.333  2.000     nano    4.18879    0.8290772
## 12 Bag     vent 36.164  3.000    micro   14.13717    2.5980292
## 13 Bag      BSW 16.206 14.924     nano 1740.42111  238.4669404
## 14 Bag      BSW  7.000  7.000     nano  179.59438   28.2640658
## 15 Bag     vent 10.069  5.000     nano   65.44985   10.9544849

Volume is reported as um^3

# Volume by experiment type
biov_calc %>% 
  group_by(EXP) %>% summarise(VOL = mean(Volume), C = mean(pgC_cell))
## # A tibble: 2 × 3
##   EXP     VOL     C
##   <chr> <dbl> <dbl>
## 1 Bag    836.  112.
## 2 IGT   3186.  400.
# Volume by euk size
biov_calc %>% 
  group_by(SizeFrac) %>% summarise(VOL = mean(Volume), C = mean(pgC_cell))
## # A tibble: 2 × 3
##   SizeFrac   VOL     C
##   <chr>    <dbl> <dbl>
## 1 micro    4678. 587. 
## 2 nano      325.  46.9
# Volume by site
biov_calc %>% 
  group_by(VENT_BSW) %>% summarise(VOL = mean(Volume), C = mean(pgC_cell))
## # A tibble: 2 × 3
##   VENT_BSW   VOL     C
##   <chr>    <dbl> <dbl>
## 1 BSW       643.  91.4
## 2 vent     2188. 276.
# head(biov_calc)
euk_vol <- mean(biov_calc$Volume);euk_vol # in um^3
## [1] 1775.759
euk_carbon <- mean(biov_calc$pgC_cell); euk_carbon # in pg C per cell
## [1] 226.9636
euk_carbon_min <- min(biov_calc$pgC_cell); euk_carbon_min
## [1] 0.8290772
euk_carbon_max <- max(biov_calc$pgC_cell); euk_carbon_max
## [1] 1120.958
# euk_carbon

Avg euk biomass pg C per individual cell == {r}euk_carbon

8.2 Protist carbon-literature

Compare with Menden-Deuer and Lessard 2000, Table 2 - using only the heterotrophic species measured. Based on Table 2, the min volume was 4745 and the maximum was 1.2 x10^7 µm^3. Carbon content was measured at pg per cell, this was 469.48-35,339 pg per cell.

Import the heterotroph species volume and carbon content to compare to my measured values.

# Hu-measured
range(biov_calc$Volume)
## [1]    4.18879 9046.50993
range(biov_calc$pgC_cell)
## [1]    0.8290772 1120.9583343
c_prev <- read.delim("input-data/md-lessard-2000.txt") # Table 2, heterotrophs only

# combine and plot
carbon_compare <- c_prev %>% 
  add_column(source = "Menden-Deuer Lessard") %>% 
  select(source, Volume = vol, pgC_cell) %>% 
  rbind(biov_calc %>% add_column(source = "MCR") %>% select(source, Volume, pgC_cell)) %>% 
    ggplot(aes(x = Volume, y = pgC_cell, fill = source)) +
      geom_point(aes(fill = source), shape = 23, color = "black", size = 3) +
      scale_y_log10() + scale_x_log10() +
      labs(title = "Compare literature to measured cell volume & C content",
           x = bquote("Volume" ~µm^-3),
           y = bquote("pg C" ~cell^-1)) +
      theme_bw() + theme(legend.title = element_blank(),
                         axis.title = element_text(size = 14),
                         axis.text = element_text(size = 14),
                         legend.text = element_text(size = 14))

carbon_compare

Upon comparison, the measured carbon content was much lower from the grazing experiments. This makes sense, as I am looking at preserved specimen and a smaller total number of cells. AND the deep-sea protist cell sizes may be smaller overall.

Find lowest estimates or protist carbon, benthic estimates, and others? How does it compare to my measurements?

euk_carbon_lit_mean <- mean(c_prev$pgC_cell)
euk_carbon_lit_min <- min(c_prev$pgC_cell)
euk_carbon_lit_max <- max(c_prev$pgC_cell)

8.3 Prokaryote carbon-literature

Below adding in biomass estimates from prokaryotes and protists.

bac_carbon_ug <- (86)*(1.00E-9) # From Derived from Morono et al. 2011 
# bac_carbon_ug
bac_carbon_ug_2 <- (173)*(1.00E-9) # Derived from McNichol et al. 2018; LOFERER-KRO ̈ ßBACHER, J. KLIMA & R. PSENNER 1998
# table_wcalcs

8.4 Generate carbon biomass table

Incorporate calculations that include biomass of population and ug C consumed. For rate measurements, only incorporate the Morono et al. 2011 biomass for prokaryotes. This way it is on the lower end and is comparable to Gorda Ridge work.

NEED TO MODIFY SO IT INCLUDES BIOMASS OF EUK AND PROK AT EACH SITE

bsw <- c("Plume", "Background")

table_wcalcs_biomass <- table_wcalcs %>% 
  add_column(euk_C_ug_Hu = (euk_carbon / (1.00E+06))) %>% # Convert to ug from pg
  add_column(euk_C_ug_lit = (euk_carbon_lit_mean / (1.00E+06))) %>% # literature
  add_column(bac_C_ug = bac_carbon_ug) %>% 
  add_column(bac_C_ug_2 = bac_carbon_ug_2) %>%
  # Grazing rate in ug C per bac per day
  mutate(RATE_ugCbac_pergrazer_perday = (RATE_hr * 24 * bac_C_ug), # Grazing rate as ug C per grazer per day
         # % of cell carbon per day
         SPEC_INGESTION_RATE = (RATE_ugCbac_pergrazer_perday / euk_C_ug_Hu),
         SPEC_INGESTION_RATE_lit = (RATE_ugCbac_pergrazer_perday / euk_C_ug_lit),
         Prok_biomass = PROK_ml * bac_carbon_ug,
         Euk_biomass_Hu = EUK_ml * euk_C_ug_Hu,
         Euk_biomass_lit = EUK_ml * euk_C_ug_lit,
         Prok_biomass_L = PROK_ml * bac_carbon_ug * 1000,
         Euk_biomass_Hu_L = EUK_ml * euk_C_ug_Hu * 1000,
         Euk_biomass_lit_L = EUK_ml * euk_C_ug_lit * 1000,
         # Repeat with SEM values
         Prok_biomass_sem = PROK_sem * bac_carbon_ug,
         Euk_biomass_Hu_sem = EUK_sem * euk_C_ug_Hu,
         Euk_biomass_lit_sem = EUK_sem * euk_C_ug_lit,
         Prok_biomass_sem_L = PROK_sem * (bac_carbon_ug* 1000),
         Euk_biomass_Hu_sem_L = EUK_sem * (euk_C_ug_Hu * 1000),
         Euk_biomass_lit_sem_L = EUK_sem * (euk_C_ug_lit * 1000)) %>% 
  type.convert(as.is = TRUE) %>%
  mutate(detected = case_when(
    RATE_min < 0 ~ "Not detected",
    TRUE ~ "Detected")) %>% 
  mutate(type = case_when(
    NAME %in% bsw ~ NAME,
    TRUE ~ paste("Vent", EXP, sep="-")
  )) %>% 
  mutate(GRAZE_RATE = case_when(
    RATE_min < 0 ~ 0,
    TRUE ~ RATE_min
  )) %>% 
  mutate(type_site = case_when(
    NAME %in% bsw ~ NAME,
    TRUE ~ "Vent"
  )) %>%
  data.frame
# View(table_wcalcs_biomass)
  • Grazing rate column == FLP per minute consumed
  • Grazing effect hr == cells per ml per hr

Also make a “bounded” table that demonstrates the ug C consumed in the context of McNichol et al.

# G = number of cells grazed during experiment duration
table_wcalcs_biomass_bounded <- table_wcalcs_biomass %>% 
  add_column(fgC_cell = 86) %>% # Add in Morono et al. 2011 value
  mutate(
    # cells_consumed_perday = (G / 1), # Rate of cells consumed * in situ prok, per day
    fgC_ml_perday = (GRAZING_EFFECT_day * fgC_cell), # Convert cell amount to fg C
    ugC_L_perday = (fgC_ml_perday * (1e-09) * 1000), # Convert to ug C per L
    lower_mcnichol = 100*(ugC_L_perday / 17.3),
    upper_mcnichol = 100*(ugC_L_perday / 321.4)
  ) %>% 
  data.frame
head(table_wcalcs_biomass_bounded)
##       SITE         NAME EXP   SAMPLE    RATE_min REPs    FLP_ml   PROK_ml
## 1 Von Damm   Background Bag   VD-BSW 0.002958889    3  4861.824  37889.62
## 2 Von Damm        Plume Bag VD-Plume 0.005274231    3 34242.354  16478.31
## 3 Von Damm         X-18 Bag   VD-X18 0.001744429    2  3148.722 111429.78
## 4 Von Damm Old Man Tree IGT   VD-OMT 0.014510943    2  5352.828  71147.40
## 5 Von Damm   Ravelin #2 Bag  VD-Rav2 0.003470217    3 51890.942  71147.40
## 6 Von Damm   Ravelin #2 IGT  VD-Rav2 0.015395240    2  5352.828  71147.40
##   PROK_sem    EUK_ml   EUK_sem   RATE_hr  RATE_day    FLP_L CLEARANCE_RATE_ml
## 1 8608.427  91.83773  21.86613 0.1775333  4.260800  4861824      3.651579e-05
## 2 2623.935 157.77468  67.09859 0.3164538  7.594892 34242354      9.241591e-06
## 3 2973.793 314.87222 104.95741 0.1046657  2.511977  3148722      3.324070e-05
## 4       NA 472.30833 122.45031 0.8706566 20.895758  5352828      1.626536e-04
## 5       NA 409.33389  73.47019 0.2082130  4.997113 51890942      4.012512e-06
## 6       NA 620.99799 123.17702 0.9237144 22.169145  5352828      1.725657e-04
##   CLEARANCE_RATE_nL SPEC_GRAZE_RATE_hr GRAZE_RATE_DAY GRAZING_EFFECT_hr
## 1      3.651579e-11          1.3835695      33.205668         127.06388
## 2      9.241591e-12          0.1522858       3.654860          24.02685
## 3      3.324070e-11          3.7040034      88.896081        1166.28778
## 4      1.626536e-10         11.5723785     277.737084        5465.73080
## 5      4.012512e-12          0.2854798       6.851515         116.85656
## 6      1.725657e-10         12.2775993     294.662382        7624.36451
##   GRAZING_EFFECT_hr_min GRAZING_EFFECT_hr_max GRAZING_EFFECT_day
## 1              96.81058             157.31719          3049.5332
## 2              13.80868              34.24501           576.6444
## 3             777.52519            1555.05037         27990.9067
## 4            4048.68948            6882.77212        131177.5392
## 5              95.88230             137.83081          2804.5574
## 6            6112.04639            9136.68264        182984.7483
##   BAC_TURNOVER_PERC BAC_TURNOVER_PERC_min BAC_TURNOVER_PERC_max  euk_C_ug_Hu
## 1          8.048465             10.414647              6.558411 0.0002269636
## 2          3.499414              4.162182              3.018725 0.0002269636
## 3         25.119772             25.808540             24.466811 0.0002269636
## 4        184.374334                    NA                    NA 0.0002269636
## 5          3.941897                    NA                    NA 0.0002269636
## 6        257.191065                    NA                    NA 0.0002269636
##   euk_C_ug_lit bac_C_ug bac_C_ug_2 RATE_ugCbac_pergrazer_perday
## 1   0.01147298  8.6e-08   1.73e-07                 3.664288e-07
## 2   0.01147298  8.6e-08   1.73e-07                 6.531607e-07
## 3   0.01147298  8.6e-08   1.73e-07                 2.160300e-07
## 4   0.01147298  8.6e-08   1.73e-07                 1.797035e-06
## 5   0.01147298  8.6e-08   1.73e-07                 4.297517e-07
## 6   0.01147298  8.6e-08   1.73e-07                 1.906547e-06
##   SPEC_INGESTION_RATE SPEC_INGESTION_RATE_lit Prok_biomass Euk_biomass_Hu
## 1         0.001614483            3.193840e-05  0.003258508     0.02084382
## 2         0.002877822            5.693033e-05  0.001417135     0.03580910
## 3         0.000951827            1.882946e-05  0.009582961     0.07146452
## 4         0.007917726            1.566319e-04  0.006118676     0.10719678
## 5         0.001893483            3.745771e-05  0.006118676     0.09290388
## 6         0.008400232            1.661770e-04  0.006118676     0.14094392
##   Euk_biomass_lit Prok_biomass_L Euk_biomass_Hu_L Euk_biomass_lit_L
## 1        1.053653       3.258508         20.84382          1053.653
## 2        1.810146       1.417135         35.80910          1810.146
## 3        3.612524       9.582961         71.46452          3612.524
## 4        5.418786       6.118676        107.19678          5418.786
## 5        4.696281       6.118676         92.90388          4696.281
## 6        7.124700       6.118676        140.94392          7124.700
##   Prok_biomass_sem Euk_biomass_Hu_sem Euk_biomass_lit_sem Prok_biomass_sem_L
## 1     0.0007403247        0.004962814           0.2508697          0.7403247
## 2     0.0002256584        0.015228935           0.7698210          0.2256584
## 3     0.0002557462        0.023821507           1.2041746          0.2557462
## 4               NA        0.027791758           1.4048704                 NA
## 5               NA        0.016675055           0.8429222                 NA
## 6               NA        0.027956696           1.4132080                 NA
##   Euk_biomass_Hu_sem_L Euk_biomass_lit_sem_L detected       type  GRAZE_RATE
## 1             4.962814              250.8697 Detected Background 0.002958889
## 2            15.228935              769.8210 Detected      Plume 0.005274231
## 3            23.821507             1204.1746 Detected   Vent-Bag 0.001744429
## 4            27.791758             1404.8704 Detected   Vent-IGT 0.014510943
## 5            16.675055              842.9222 Detected   Vent-Bag 0.003470217
## 6            27.956696             1413.2080 Detected   Vent-IGT 0.015395240
##    type_site fgC_cell fgC_ml_perday ugC_L_perday lower_mcnichol upper_mcnichol
## 1 Background       86     262259.86   0.26225986      1.5159529     0.08159921
## 2      Plume       86      49591.42   0.04959142      0.2866556     0.01542981
## 3       Vent       86    2407217.98   2.40721798     13.9145548     0.74897884
## 4       Vent       86   11281268.37  11.28126837     65.2096438     3.51003994
## 5       Vent       86     241191.93   0.24119193      1.3941730     0.07504416
## 6       Vent       86   15736688.36  15.73668836     90.9635165     4.89629383
# View(table_wcalcs_biomass_bounded)
# write_delim(table_wcalcs_biomass_bounded, file = "table-wcalc.txt", delim = "\t")
# ?write_delim()

8.5 Compilation of plots

Format and factor values to plot. * Grazing rate column == FLP per minute consumed * Grazing effect hr == cells per ml per hr * Specific ingestion rate == % of cell carbon per day, estimated with literature and measured carbon values

biomass_rate_plot <- table_wcalcs_biomass_bounded %>% 
  select(SITE, NAME, EXP, SAMPLE, type, 
        PROK_ml, EUK_ml, PROK_sem, EUK_sem,
        Prok_biomass_L, Euk_biomass_Hu_L, Euk_biomass_lit_L,
        Prok_biomass_sem_L, Euk_biomass_Hu_sem_L, Euk_biomass_lit_sem_L,
        GRAZING_EFFECT_hr, GRAZING_EFFECT_hr_min, GRAZING_EFFECT_hr_max, BAC_TURNOVER_PERC, BAC_TURNOVER_PERC_min, BAC_TURNOVER_PERC_max,
        GRAZE_RATE, RATE_ugCbac_pergrazer_perday,
        SPEC_INGESTION_RATE, SPEC_INGESTION_RATE_lit) %>% 
  pivot_longer(cols = c(PROK_ml, EUK_ml,
                        Prok_biomass_L, Euk_biomass_Hu_L, Euk_biomass_lit_L,
                        GRAZING_EFFECT_hr, BAC_TURNOVER_PERC,
                        GRAZE_RATE, RATE_ugCbac_pergrazer_perday,
                        SPEC_INGESTION_RATE, SPEC_INGESTION_RATE_lit),
               names_to = "Variable", values_to = "Value") %>% 
  pivot_longer(cols = c(PROK_sem, EUK_sem, Prok_biomass_sem_L, Euk_biomass_Hu_sem_L, Euk_biomass_lit_sem_L,
                        BAC_TURNOVER_PERC_max, BAC_TURNOVER_PERC_min, GRAZING_EFFECT_hr_max, GRAZING_EFFECT_hr_min, GRAZING_EFFECT_hr_max),
               names_to = "SEM_variable", values_to = "SEM") %>% 
  data.frame

biomass_rate_plot$NAME_ORDER <- factor(biomass_rate_plot$NAME, levels = c("Background","Plume", "Quakeplume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole", "Hot Chimlet #1", "Shrimp Gulley", "South of Hot Chimlet", "South of LungSnack", "Arrow Loop", "Bartizan", "Ravelin #1"))

biomass_rate_plot$VARIABLE_ORDER <- factor(biomass_rate_plot$Variable, 
                                           levels = c("PROK_ml", "EUK_ml",
                                                      "Prok_biomass_L", "Euk_biomass_Hu_L", "Euk_biomass_lit_L",
                                                      "GRAZING_EFFECT_hr", "BAC_TURNOVER_PERC",
                                                      "GRAZE_RATE", "RATE_ugCbac_pergrazer_perday",
                                                      "SPEC_INGESTION_RATE", "SPEC_INGESTION_RATE_lit"),
                                           labels = c("Prokaryote~cells~mL^{-1}", "Eukaryote~cells~mL^{-1}",
                                                      "Prokaryote~µg~C~L^{-1}", "Measured~eukaryote~µg~C~L^{-1}", "Literature-based~eukaryote~µg~C~L^{-1}",
                                                      "Cells~mL^{-1}~hr^{-1}", "Prokaryote~turnover~'%'~d^{-1}",
                                                      "FLP~consumed~min^{-1}", "µg~C~grazer^{-1}~day^{-1}",
                                                      "Cell~carbon~%~day^{-1}", "Cell~carbon~%~day^{-1}~(lit)"))

# View(biomass_rate_plot)

Function to generate consistent plots for Mid-Cayman Rise cell concentration and biomass data.

conc_rate_plot_mcr <- function(df, var, sem){
  df %>% 
    filter(Variable == var) %>%
    filter(SEM_variable == sem) %>% 
    ggplot(aes(y = Value, x = NAME_ORDER, shape = EXP, fill = SITE)) +
    geom_errorbar(aes(ymax = (Value + SEM), ymin = (Value - SEM)), 
                  width = 0.2, position = position_dodge(width = 0.4)) +
    geom_point(stat = "identity", aes(shape = EXP, fill = SITE),
               color = "black", size = 3, position = position_dodge(width = 0.4)) +
    scale_shape_manual(values = c(21, 23)) +
    scale_fill_manual(values = c("#de2d26", "#1c9099")) +
    facet_wrap(VARIABLE_ORDER ~ ., scales = "free", 
               strip.position = c("left"), labeller = label_parsed) +
    scale_y_log10() +
    # scale_y_log10(labels = function(x) format(x, scientific = TRUE)) +
    theme_minimal() +
    theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
          panel.background = element_blank(), 
          axis.line = element_line(colour = "black"), 
          axis.text.x = element_text(color="black", size = 11, 
                                     angle = 45, hjust = 1, vjust = 1), 
          axis.text.y = element_text(color="black", size = 11),
          axis.title =element_text(color="black", size = 14),
          axis.ticks = element_line(),
          legend.title = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(color="black", size = 11),
          strip.text.x = element_blank())+
    guides(fill = guide_legend(override.aes = list(shape = c(21))),
           shape = guide_legend(override.aes = list(fill = "black"))) +
    labs(x = "", y = "")
}
# svg("cells-per-ml-yaxies.svg", h = 8, w = 13)
plot_grid(conc_rate_plot_mcr(biomass_rate_plot, "PROK_ml", "PROK_sem"),
          conc_rate_plot_mcr(biomass_rate_plot, "EUK_ml", "EUK_sem"),
          conc_rate_plot_mcr(biomass_rate_plot, "Prok_biomass_L", "Prok_biomass_sem_L"),
          conc_rate_plot_mcr(biomass_rate_plot, "Euk_biomass_Hu_L", "Euk_biomass_Hu_sem_L"))

# dev.off()
# svg("MCR-only-cellcounts-grazingrate.svg", h=12, w=7)
plot_grid(conc_rate_plot_mcr(biomass_rate_plot, "PROK_ml", "PROK_sem"),
          conc_rate_plot_mcr(biomass_rate_plot, "EUK_ml", "EUK_sem"),
          conc_rate_plot_mcr(biomass_rate_plot, "GRAZING_EFFECT_hr", "Prok_biomass_sem_L"), ncol = 1)
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning in self$trans$transform(x): NaNs produced
## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

# dev.off()
subset <- c("GRAZING_EFFECT_hr", "RATE_ugCbac_pergrazer_perday")
biomass_rate_plot %>%
  filter(Variable %in% subset) %>% 
  group_by(EXP, type, Variable) %>% 
  summarise(MEAN = mean(Value))
## `summarise()` has grouped output by 'EXP', 'type'. You can override using the `.groups` argument.
## # A tibble: 8 × 4
## # Groups:   EXP, type [4]
##   EXP   type       Variable                        MEAN
##   <chr> <chr>      <chr>                          <dbl>
## 1 Bag   Background GRAZING_EFFECT_hr            1.27e+2
## 2 Bag   Background RATE_ugCbac_pergrazer_perday 3.66e-7
## 3 Bag   Plume      GRAZING_EFFECT_hr            3.41e+1
## 4 Bag   Plume      RATE_ugCbac_pergrazer_perday 6.59e-7
## 5 Bag   Vent-Bag   GRAZING_EFFECT_hr            1.21e+3
## 6 Bag   Vent-Bag   RATE_ugCbac_pergrazer_perday 4.31e-7
## 7 IGT   Vent-IGT   GRAZING_EFFECT_hr            8.58e+3
## 8 IGT   Vent-IGT   RATE_ugCbac_pergrazer_perday 1.50e-6
# 8.577791e+03  - 1.214981e+03  
# scale_y_continuous(labels = function(x) format(x, scientific = TRUE)) +
grazing_rate <- biomass_rate_plot %>% 
  type.convert(as.is = TRUE) %>%
    filter(Variable == "GRAZING_EFFECT_hr") %>%
    filter(SEM_variable == "GRAZING_EFFECT_hr_min" | SEM_variable == "GRAZING_EFFECT_hr_max") %>% 
    pivot_wider(names_from = SEM_variable, values_from = SEM) %>% 
    ggplot(aes(y = Value, x = NAME_ORDER, shape = EXP, fill = SITE)) +
    geom_errorbar(aes(ymax = (GRAZING_EFFECT_hr_max), ymin = (GRAZING_EFFECT_hr_min)), 
                  width = 0.2, position = position_dodge(width = 0.4)) +
    geom_point(stat = "identity", aes(shape = EXP, fill = SITE),
               color = "black", size = 3, position = position_dodge(width = 0.4)) +
    scale_shape_manual(values = c(21, 24)) +
    scale_fill_manual(values = c("#de2d26", "#1c9099")) +
    facet_wrap(VARIABLE_ORDER ~ ., scales = "free", 
               strip.position = c("left"), labeller = label_parsed) +
    theme_minimal() +
    theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
          panel.background = element_blank(), 
          axis.line = element_line(colour = "black"), 
          axis.text.x = element_text(color="black", size = 11, 
                                     angle = 45, hjust = 1, vjust = 1), 
          axis.text.y = element_text(color="black", size = 11),
          axis.title =element_text(color="black", size = 11),
          axis.ticks = element_line(),
          legend.title = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(color="black", size = 11),
          strip.text.x = element_blank())+
    guides(fill = guide_legend(override.aes = list(shape = c(21))),
           shape = guide_legend(override.aes = list(fill = "black"))) +
    labs(x = "", y = "")
# class(grazing_rate)
# head(biomass_rate_plot)
# svg("plot-grazing-withscale.svg", h = 7, w = 6)
# grazing_rate
plot_grid(grazing_rate + scale_y_continuous(limits = c(0,2000)),
          grazing_rate,
          ncol = 1, rel_heights = c(0.6,1))
## Warning: Removed 4 rows containing missing values (geom_point).

# dev.off()
# scale_y_continuous(labels = function(x) format(x, scientific = TRUE)) +
prok_turnover <- biomass_rate_plot %>% 
    filter(Variable == "BAC_TURNOVER_PERC") %>%
    filter(SEM_variable == "BAC_TURNOVER_PERC_min" | SEM_variable == "BAC_TURNOVER_PERC_max") %>% 
    pivot_wider(names_from = SEM_variable, values_from = SEM) %>% 
    ggplot(aes(y = Value, x = NAME_ORDER, fill = SITE)) +
    geom_bar(stat = "identity", aes(fill = SITE),
               color = "black", width = 2, position = position_dodge2(preserve = "single")) +
    geom_errorbar(aes(ymax = (BAC_TURNOVER_PERC_max), ymin = (BAC_TURNOVER_PERC_min)), 
                  width = 0.3, position = position_dodge2(preserve = "single")) +
    scale_fill_manual(values = c("#de2d26", "#1c9099")) +
    facet_wrap(VARIABLE_ORDER ~ ., scales = "free", 
               strip.position = c("left"), labeller = label_parsed) +
    theme_minimal() +
    theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
          panel.background = element_blank(), 
          axis.line = element_line(colour = "black"), 
          axis.text.x = element_text(color="black", size = 11, 
                                     angle = 45, hjust = 1, vjust = 1), 
          axis.text.y = element_text(color="black", size = 11),
          axis.title =element_text(color="black", size = 11),
          axis.ticks = element_line(),
          legend.title = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(color="black", size = 11),
          strip.text.x = element_blank())+
    guides(fill = guide_legend(override.aes = list(shape = c(21))),
           shape = guide_legend(override.aes = list(fill = "black"))) +
    labs(x = "", y = "")
# View(biomass_rate_plot)
# tmp <- (biomass_rate_plot %>% filter(Variable == "RATE_ugCbac_pergrazer_perday" & Value > 0))
# range(tmp$Value)
# View(tmp)
# unique(biomass_rate_plot$Variable)
# svg("prok-turnover.svg", h = 8, w = 6)
plot_grid(prok_turnover,
          prok_turnover + scale_y_continuous(limits = c(0,100)),
          ncol = 1)
## Warning: Removed 3 rows containing missing values (geom_bar).

# dev.off()

8.6 Compare MCR data with Gorda Ridge - think about environmental parameters

Import Gorda Ridge data

gr <- read.delim("../../GordaRidgeCruise_2019/protist-gordaridge-2021/Grazing-at-GordaRidge-SKH-2021/data-input/Grazing-calc-wCarbon-results.txt")
# env_gr <- read.delim("../../GordaRidgeCruise_2019/protist-gordaridge-2021/Grazing-at-GordaRidge-SKH-2021/data-input/GR-environ-SAMPLE.txt")

temps <- read.delim("/Users/sarahhu/Desktop/Projects/GordaRidgeCruise_2019/protist-gordaridge-2019-analysis/temperature-allvents.txt")
mcr_graze <- read.delim("table-wcalc.txt")
# head(mcr_graze)
all_vents <- mcr_graze %>%
  type.convert(as.is = TRUE) %>%
  select(SITE, NAME, SAMPLE, EXP, PROK_ml, EUK_ml, GRAZING_EFFECT_hr, GRAZING_EFFECT_hr_min, GRAZING_EFFECT_hr_max, BAC_TURNOVER_PERC, ugC_L_perday) %>% 
  rbind(gr %>%
          add_column(SITE = "Gorda Ridge") %>% 
          add_column(EUK_ml = NA) %>% 
          separate(SAMPLE, c("SAMPLE", "NAME"), sep = "-") %>% 
          select(SITE, NAME, SAMPLE, EXP = Bottle, PROK_ml = prok_avg, EUK_ml, GRAZING_EFFECT_hr = GrazingRate_hr, GRAZING_EFFECT_hr_min = GrazingRate_hr_min, GRAZING_EFFECT_hr_max = GrazingRate_hr_max, BAC_TURNOVER_PERC = Prok_turnover, ugC_L_perday = ugC_L_perday_morono)) %>% 
  left_join(temps) %>% 
  mutate(SAMPLE_TYPE = case_when(
    grepl("BSW", NAME) ~ "Background",
    grepl("Near vent BW", NAME) ~ "Background",
    grepl("Background", NAME) ~ "Background", 
    grepl("Plume", NAME) ~ "Background",
    TRUE ~ "Vent"
  ))
## Joining, by = c("NAME", "SAMPLE")
# View(all_vents)
# all_vents
# write_delim(all_vents, file = "grazing-cellcounts-GR_MCR.txt", delim = "\t")

Plot grazing rates

allrates <- all_vents %>% 
  select(SITE, NAME, SAMPLE, EXP, SAMPLE_TYPE, starts_with("GRAZING_EFFECT_")) %>% 
  distinct() %>% 
    ggplot(aes(y = GRAZING_EFFECT_hr, x = NAME, fill = SITE, shape = SAMPLE_TYPE)) +
    geom_errorbar(aes(ymax = (GRAZING_EFFECT_hr_max), ymin = (GRAZING_EFFECT_hr_min)), 
                  width = 0.2, position = position_dodge(width = 0.4)) +
    geom_point(stat = "identity", aes(fill = SITE, shape = SAMPLE_TYPE),
               color = "black", size = 3, position = position_dodge(width = 0.4)) +
    scale_shape_manual(values = c(21, 24)) +
    scale_fill_manual(values = c("#de2d26", "#1c9099", "#addd8e")) +
    facet_grid(. ~ SAMPLE_TYPE, scales = "free", space = "free") +
    theme_minimal() +
    theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
          panel.background = element_blank(), 
          axis.line = element_line(colour = "black"), 
          axis.text.x = element_text(color="black", size = 11, 
                                     angle = 45, hjust = 1, vjust = 1), 
          axis.text.y = element_text(color="black", size = 11),
          axis.title =element_text(color="black", size = 11),
          axis.ticks = element_line(),
          legend.title = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(color="black", size = 11),
          strip.text.x = element_blank())+
    guides(fill = guide_legend(override.aes = list(shape = c(21))),
           shape = guide_legend(override.aes = list(fill = "black"))) +
    labs(x = "", y = bquote("cells"~ml^-1~hr^-1))
# svg("compare-all-rates.svg", h = 10, w = 7)
plot_grid(allrates,
          allrates + scale_y_continuous(limits = c(0,1000)),
          allrates + scale_y_log10(), ncol = 1)
## Warning: Removed 7 rows containing missing values (geom_point).
## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

# dev.off()
# svg("compare-all-rates-color.svg", h = 4, w = 9)
# allrates + scale_y_log10()
# dev.off()

Repeat grazing rate plot, but removed undetectable

# unique(all_vents$GRAZING_EFFECT_hr)
allrates_nonzero <- all_vents %>% 
  filter(GRAZING_EFFECT_hr > 0) %>% 
  select(SITE, NAME, SAMPLE, EXP, SAMPLE_TYPE, starts_with("GRAZING_EFFECT_")) %>% 
  distinct() %>% 
    ggplot(aes(y = GRAZING_EFFECT_hr, x = NAME, fill = SITE, shape = SAMPLE_TYPE)) +
    geom_errorbar(aes(ymax = (GRAZING_EFFECT_hr_max), ymin = (GRAZING_EFFECT_hr_min)), 
                  width = 0.2, position = position_dodge(width = 0.4)) +
    geom_point(stat = "identity", aes(fill = SITE, shape = SAMPLE_TYPE),
               color = "black", size = 3, position = position_dodge(width = 0.4)) +
    scale_shape_manual(values = c(21, 24)) +
    scale_fill_manual(values = c("#de2d26", "#1c9099", "#addd8e")) +
    facet_grid(. ~ SAMPLE_TYPE, scales = "free", space = "free") +
    theme_minimal() +
    theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
          panel.background = element_blank(), 
          axis.line = element_line(colour = "black"), 
          axis.text.x = element_text(color="black", size = 11, 
                                     angle = 45, hjust = 1, vjust = 1), 
          axis.text.y = element_text(color="black", size = 11),
          axis.title =element_text(color="black", size = 11),
          axis.ticks = element_line(),
          legend.title = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(color="black", size = 11),
          strip.text.x = element_blank())+
    guides(fill = guide_legend(override.aes = list(shape = c(21))),
           shape = guide_legend(override.aes = list(fill = "black"))) +
    labs(x = "", y = bquote("cells"~ml^-1~hr^-1))
# svg("compare-all-rates-color-nonZero.svg", h = 4, w = 7)
allrates_nonzero + scale_y_log10()

# dev.off()
# svg("plot-gcellconc-grazing-log.svg", h = 10, w = 6)
# grazing_rate + scale_y_log10(limits = c(1e-1,1e5), breaks=10^(0:7))
plot_grid(conc_rate_plot_mcr(biomass_rate_plot, "PROK_ml", "PROK_sem"),
          conc_rate_plot_mcr(biomass_rate_plot, "EUK_ml", "EUK_sem"),
allrates_nonzero + scale_y_log10(), ncol = 1)

# dev.off()

8.7 Run regression analysis with MCR and GR data

metadata <- read.delim("input-data/flp-exp-metadata-compiled.txt")
# head(metadata)
# temps
# all_vents

Linear regression with both Gorda Ridge and MCR data

library(broom)
# ?pivot_longer
# e IGT
# head(all_vents)
regression_input <- all_vents %>% 
  filter(!(SAMPLE == "Piccard-Shrimpocalypse" & EXP == "IGT")) %>% 
  filter(!is.na(Highest.Temp)) %>% 
    select(SITE, NAME, SAMPLE, SAMPLE_TYPE, EXP, PROK_ml, EUK_ml, GRAZING_EFFECT_hr, BAC_TURNOVER_PERC, ugC_L_perday, TEMP = Highest.Temp) %>% 
  mutate(PROK_EUK_RATIO = (PROK_ml/EUK_ml)) %>% 
  pivot_longer(cols = c(GRAZING_EFFECT_hr, BAC_TURNOVER_PERC, ugC_L_perday), names_to = "RATE", values_to = "RATE_VALUE") %>% 
  pivot_longer(cols = c(PROK_ml, EUK_ml, PROK_EUK_RATIO, TEMP), names_to = "PARAMS", values_to = "PARAMS_VALUE") %>% 
  data.frame

regression_tmp <- regression_input %>% 
  # Set up the linear regression
  group_by(RATE, PARAMS) %>% 
  nest(-RATE, -PARAMS) %>% 
  mutate(lm_fit = map(data, ~lm(RATE_VALUE ~ PARAMS_VALUE, data = .)), 
    tidied = map(lm_fit, tidy)) %>% 
  unnest(tidied) %>% 
  select(RATE, PARAMS, 
    term, estimate) %>% 
  pivot_wider(names_from = term, values_from = estimate) %>% 
  select(everything(), SLOPE = PARAMS_VALUE) %>%
  data.frame
## Warning: All elements of `...` must be named.
## Did you want `data = c(SITE, NAME, SAMPLE, SAMPLE_TYPE, EXP, RATE_VALUE, PARAMS_VALUE)`?
regression_results <- regression_input %>% 
  group_by(RATE, PARAMS) %>% 
  nest(-RATE, -PARAMS) %>% 
  mutate(lm_fit = map(data, ~lm(RATE_VALUE ~ PARAMS_VALUE, data = .)),
    glanced = map(lm_fit, glance)) %>% 
  unnest(glanced) %>% 
  select(RATE, PARAMS, r.squared, adj.r.squared) %>% 
  right_join(regression_tmp) %>%
  right_join(regression_input) %>%
  data.frame
## Warning: All elements of `...` must be named.
## Did you want `data = c(SITE, NAME, SAMPLE, SAMPLE_TYPE, EXP, RATE_VALUE, PARAMS_VALUE)`?
## Joining, by = c("RATE", "PARAMS")
## Joining, by = c("RATE", "PARAMS")
# head(regression_results)

Plot regression results

regression_results %>% 
  # filter(RATE == "GRAZING_EFFECT_hr") %>%
  # filter(PARAMS == "TEMP") %>% 
  ggplot(aes(x = PARAMS_VALUE, y = RATE_VALUE, shape = SAMPLE_TYPE, fill = SITE)) +
    geom_abline(aes(slope = SLOPE, intercept = `X.Intercept.`), color = "black", linetype = "dashed", size = 0.5) +
    geom_point(color = "black", aes(shape = SAMPLE_TYPE, fill = SITE)) +
    scale_shape_manual(values = c(21, 24)) +
    scale_fill_manual(values = c("#476AA7","#7299CE", "#A2937A")) +
    facet_wrap(PARAMS ~ RATE + round(r.squared, 3), scales = "free", ncol = 5,
             strip.position = "bottom", labeller = label_parsed) +
  theme_bw() +
  theme(
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text = element_text(color = "black", size = 10),
    axis.title = element_text(color = "black", size = 10),
    legend.title = element_blank()) # +
## Warning: Removed 24 rows containing missing values (geom_point).

  # labs(y = bquote("Cells "~mL^-1 ~hr^-1), x = "Temperature (C)")

NEXT STEP = get more biogeochem data?

9 Parameter comparions to plot

head(regression_results)
##                RATE  PARAMS r.squared adj.r.squared X.Intercept.      SLOPE
## 1 GRAZING_EFFECT_hr PROK_ml 0.3260524     0.2864084    -488.8238 0.02782789
## 2 GRAZING_EFFECT_hr PROK_ml 0.3260524     0.2864084    -488.8238 0.02782789
## 3 GRAZING_EFFECT_hr PROK_ml 0.3260524     0.2864084    -488.8238 0.02782789
## 4 GRAZING_EFFECT_hr PROK_ml 0.3260524     0.2864084    -488.8238 0.02782789
## 5 GRAZING_EFFECT_hr PROK_ml 0.3260524     0.2864084    -488.8238 0.02782789
## 6 GRAZING_EFFECT_hr PROK_ml 0.3260524     0.2864084    -488.8238 0.02782789
##       SITE         NAME   SAMPLE SAMPLE_TYPE EXP RATE_VALUE PARAMS_VALUE
## 1 Von Damm   Background   VD-BSW  Background Bag  127.06388     37889.62
## 2 Von Damm   Background   VD-BSW  Background Bag  127.06388     37889.62
## 3 Von Damm   Background   VD-BSW  Background Bag  127.06388     37889.62
## 4 Von Damm        Plume VD-Plume  Background Bag   24.02685     16478.31
## 5 Von Damm         X-18   VD-X18        Vent Bag 1166.28778    111429.78
## 6 Von Damm Old Man Tree   VD-OMT        Vent IGT 5465.73080     71147.40
unique(regression_results$PARAMS)
## [1] "PROK_ml"        "EUK_ml"         "PROK_EUK_RATIO" "TEMP"

Plot temperature by euk:prok ratio? individual? other things?

10 Prediction of MCR data with parameters

From Vaque et al. 1994 paper equation (1):

log(GT) = -3.21 + 0.99 log (HNF) + 0.028 (T) + 0.55 log(BAC)

GT = grazing rate (cells ml-1 hr-1) HNF = heterotrophic nanoflagellates (cells ml-1) T = temperature (C) BAC = prokaryote abundance (cells ml-1)

Start with non-zero grazing rate values, etc. Where… GT == Grazing effect hr HNF == Euk_ml BAC == PROK_ml T == highest temp

We only have these for MCR data.. can include GR for the “predicting euk_ml”

10.1 Predicted parameters

# head(all_vents)
allvent_wpredicted <- all_vents %>% 
  select(SITE, NAME, SAMPLE, EXP, PROK_ml, EUK_ml, GRAZING_EFFECT_hr, Highest.Temp) %>% 
  group_by(SITE, NAME, SAMPLE, EXP, PROK_ml, EUK_ml, GRAZING_EFFECT_hr) %>% 
  summarise(TEMP = mean(Highest.Temp)) %>% 
  distinct() %>% 
  mutate(PREDICTED_TEMP = (((log(GRAZING_EFFECT_hr)) + 3.21 - (0.99*log(EUK_ml)) - (0.55*log(PROK_ml))) / 0.028),
         PREDICTED_GRAZING_EFFECT_hr = (10^(-3.21 + (0.99*log(EUK_ml)) + (0.028*TEMP) + (0.55*(log(PROK_ml))))),
         PREDICTED_EUK = (10^(((log(GRAZING_EFFECT_hr)) + 3.21 - (0.55*log(PROK_ml)) - (0.028*TEMP))/0.99)),
         PREDICTED_PROK = (10^(((log(GRAZING_EFFECT_hr)) + 3.21 - (0.99*log(EUK_ml)) - (0.028*TEMP))/0.55)))
## `summarise()` has grouped output by 'SITE', 'NAME', 'SAMPLE', 'EXP', 'PROK_ml', 'EUK_ml'. You can override using the `.groups` argument.
# View(allvent_wpredicted)
sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS/LAPACK: /Users/sarahhu/anaconda3/envs/r_4.1/lib/libopenblasp-r0.3.15.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] gt_0.3.1        broom_0.7.9     cowplot_1.1.1   forcats_0.5.1  
##  [5] stringr_1.4.0   dplyr_1.0.7     purrr_0.3.4     readr_2.0.0    
##  [9] tidyr_1.1.3     tibble_3.1.3    ggplot2_3.3.5   tidyverse_1.3.1
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.7         lubridate_1.7.10   assertthat_0.2.1   digest_0.6.27     
##  [5] utf8_1.2.2         R6_2.5.0           cellranger_1.1.0   backports_1.2.1   
##  [9] reprex_2.0.1       evaluate_0.14      httr_1.4.2         highr_0.9         
## [13] pillar_1.6.2       rlang_0.4.11       readxl_1.3.1       rstudioapi_0.13   
## [17] jquerylib_0.1.4    checkmate_2.0.0    rmarkdown_2.9      labeling_0.4.2    
## [21] munsell_0.5.0      compiler_4.1.0     modelr_0.1.8       xfun_0.24         
## [25] pkgconfig_2.0.3    htmltools_0.5.2    tidyselect_1.1.1   fansi_0.5.0       
## [29] crayon_1.4.1       tzdb_0.1.2         dbplyr_2.1.1       withr_2.4.2       
## [33] grid_4.1.0         jsonlite_1.7.2     gtable_0.3.0       lifecycle_1.0.0   
## [37] DBI_1.1.1          magrittr_2.0.1     scales_1.1.1       cli_3.0.1         
## [41] stringi_1.7.4      farver_2.1.0       fs_1.5.0           xml2_1.3.2        
## [45] bslib_0.3.0        ellipsis_0.3.2     generics_0.1.0     vctrs_0.3.8       
## [49] RColorBrewer_1.1-2 tools_4.1.0        glue_1.4.2         hms_1.1.0         
## [53] fastmap_1.1.0      yaml_2.2.1         colorspace_2.0-2   rvest_1.0.1       
## [57] knitr_1.33         haven_2.4.3        sass_0.4.0